Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 18 January 2013 



(MN WTe^ style file v2.2) 



Edge modes in self-gravitating disc-planet interactions 



o 



Min-Kai Lin^ and John C. B. Papaloizou^ f 

^ Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, 
Wilberforce Road, Cambridge, CBS OWA, UK 



18 January 2013 



ABSTRACT 

We study the stability of gaps opened by a giant planet in a self-gravitating proto- 
planetary disc. We find a linear instability associated with both the self-gravity of the 
disc and local vortensity maxima which coincide with gap edges. For our models, these 
edge modes develop and extend to twice the orbital radius of a Saturn mass planet 
in discs with total masses Md > 0.06M* where is the central stellar mass, corre- 
sponding to a Toomre Q < 1.5 at twice the planet's orbital radius. The disc models, 
although massive, are such that they are stable in the absence of the planet. Unlike 
the previously studied local vortex forming instabilities associated with gap edges in 
weakly or non-self-gravitating discs with low viscosity, the edge modes we consider are 
global and exist only in sufficiently massive discs, but for the typical viscosity values 
adopted for protoplanetary discs. 

It is shown through analytic modelling and linear calculations that edge modes 
may be interpreted as a localised disturbance associated with a gap edge inducing 
activity in the extended disc, through the launching of density waves excited through 
gravitational potential perturbation at Lindblad resonances. We also perform hydro- 
dynamic simulations in order to investigate the evolution of edge modes in the linear 
and nonlinear regimes in disc-planet systems. The form and growth rates of developing 
unstable modes are found to be consistent with linear theory. Their dependence on 
viscosity and gravitational softening is also explored. 

We also performed a first study of the effect of edge modes on disc-planet torques 
and the orbital migration of the planet. We found that if edge modes develop, then 
the average torque on the planet becomes more positive with increasing disc mass. In 
simulations where the planet was allowed to migrate, although a fast type III migration 
could be seen that was similar to that seen in non self-gravitating discs, we found that 
it was possible for the planet to interact gravitationally with the spiral arms associated 
with an edge mode and that this could result in the planet being scattered outwards. 
Thus orbital migration is likely to be complex and non monotonic in massive discs of 
the type we consider. 
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1 INTRODUCTION 

Understanding the interaction between gaseous protoplane- 
tary discs and embedded planets is important for planet for- 
mation theory. Such interaction may lead to inward orbital 
migration and account for at least so me of the observed class 
of exo-planets called 'hot Jupiters avor k Q uelodfTooil ) 
which orbit close to their central stars. Analytical studies of 
disc-planet intera ction began well before the d iscovery of ex- 
trasolar planets (jColdreich TremaindflOTOl l. It is known 
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that giant planets with masses comparable to or exceed- 
ing t hat of Saturn can open gaps in standard model discs 
( Papaloizou k LinlfTosl ) subsequent to which they may mi- 
grate inwards. 

The stability of protoplanetary discs with gaps opened 
by int eraction with a planet has also been the subject of 
study (|Koller et al.ll2003l : iLi et~al]l2005l : Ide Val-Borro et al.1 
2007), as w ell as the consequ ences of instability on planetary 
migra tion (lOu et al. 2007: L i et alllioO^ : iLin PapaloizoT] 
bold : IYu et al.ll20T^ ). These works focused on low viscos- 
ity discs where gap edges become unstable with the result 
that vortices form. Such instabilities are known to be asso- 
ciated with steep surface density gradients or narrow rings 
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l200d . I 200 ih . These works, like most disc-planet studies, ei- 
ther ignore self-gravity completely or at best consider weak 
self-gravity. We remark that partial disc gaps induced by a 
Saturn ma ss planet may be associated with rapid type II I 
migration (jMasset k Papaloizoulliooi : IPeplihski et al.ll2008l ) 
in a massive enough disc wh ich is viscous enough to be sta- 
ble against vortex formation. Lin & Papaloizou (2010) found 
that type III migration, although unsteady, could on average 
persist when the disc viscosity was small enough for vortex 
instability to occur. 

However, as type III migration occurs in massive discs 
the effects of self-gravity need to be properly considered for 
consistency. The stability of a disc gap induced by inter- 
action with a giant planet in massive, self- gravitating discs 
has not yet been studied. The stability of structured self- 
gravitating discs without planets was explored for particle 
discs bv ISellwood KahnI (Il99lh and for gaseous discs b y 
IPapaloizou Linl(| 19891 ) andlPapal oizou &: Savonii i_(ll99lh . 
Similarly, more recently iMeschiari LaughlinI (|2008l l also 
demonstrated gravitational instabilities associated with pre- 
scribed surface density profiles that model gap structure. 

In this paper, we extend the above works by studying 
the gravitational stability of gaps self-consistently opened 
by a planet. In this respect the disc models are stable if 
the planet is not introduced and there is thus no gap. We 
extend the analytic discussion of neutral modes associated 
with surface density gap edges given by ISellwood KahnI 
(|l99lh to gaseous discs. In addition we develop the physi- 
cal interpretation of the associated gap edge instabilities in 
massive discs as disturbances localised around a vortensity 
maximum that further perturb gravitationally the smooth 
parts of the disc by exciting waves at Lindblad resonances. 
The angular momentum carried away reacts back on the 
edge disturbance so as to destabilise it. 

We perform both linear calculations and nonlinear sim- 
ulations for discs with a range of masses that consistently 
identify the growth of low azimuthal mode number, m, edge 
modes as being the dominant form of instability. In addition 
we evaluate the effect of these on the disc torques acting on 
the planet. We find that these torques become unsteady and 
oscillate in time. First estimates of the effect on the plane- 
tary migration show that although fast type Ill-like inward 
migration may occur in the unstable massive discs we study, 
scattering by spiral arms may also occur leading to short pe- 
riods of significant outward migration. 

This paper is organised as follows. We present the gov- 
erning equations and basic model in In 33]we then present 
the results of an illustrative simulation of a self- gravitating 
disc with an embedded Saturn mass planet that produces a 
dip/gap in the surface density profile, and subsequently un- 
dergoes an instability in which the gap edges play an impor- 
tant role. In the absence of the planet and gap, the Toomre 
Q value is high enough for the disc model to remain stable. 

In 311 we present an analytic discussion of the linearly 
unstable modes. We derive the general wave action conserva- 
tion law for these modes and apply it to study their angular 
momentum balance. In particular we show that, when the 
equation of state is barotropic, at marginal stability the an- 
gular momentum loss through waves propagating out of the 
system is balanced by corotation torques exerted on the disc. 
These are expected to be strongest near an edge. In 35]we 



consider modes localised around a vortensity maximum near 
an edge for which the self-gravity response balances the po- 
tentially singular response at corotation, showing that such 
disturbances may be driven by wave excitation at Lindblad 
resonances. In ^SJwe go on to present linear calculations for 
various disc models which confirm the existence of edge dom- 
inated modes for low values of the azimuthal mode number. 
In addition we find instabilities for larger values of m for 
which edge effects are less dominant but these modes have 
weaker growth and do not appear in nonlinear numerical 
simulations. 

In 33 we present results from hydrodynamic simulations 
for a range of disc masses. These are all stable in the absence 
of the planet. However, they exhibit low m edge dominated 
modes once a planetary induced gap is present. The form 
and behaviour of these is found to be in accord with linear 
theory. However, in the simulations with the lower Q values, 
the gap is found to widen and deepen as the simulation pro- 
gresses. In addition the global angular momentum transport 
through the disc, measured through an effective a parame- 
ter is increased above the level that would be induced by the 
planet alone. In addition the spiral arms associated with the 
edge instabilities are shown to produce fluctuating torques 
acting on the planet. Fast inwards type HI migration may 
occur, but we also observed outwards migration due to in- 
teraction with spiral arms. This indicates that migration is 
unlikely to be a simple monotonic process in a gravitation- 
ally active disc. Finally in ^HJwe summarise our results and 
conclude. 



2 BASIC EQUATIONS AND MODEL 

We describe here the governing equations and models for 
the disc-planet system as used in analytic discussions, lin- 
ear calculations and solved in hydrodynamic simulations. 
The system we consider is a gaseous disc of mass Md or- 
biting a central star of mass M*. We adopt a cylindrical, 
non-rotating, co-ordinate system (r, cp, z) centred on the star 
where z is the vertical co-ordinate increasing in the direc- 
tion normal to the disc mid-plane. It is convenient to adopt 
a system of equations in three dimensions to begin the dis- 
cussion of angular momentum conservation in section 14.11 
below. We begin by listing these and then go on to explain 
how we adapt them to obtain the system governing a two 
dimensional razor thin disc that we solve in hydrodynamic 
simulations. They are the continuity equation 



|+V.(,u) = 0, 

and the equation of motion 

p + u ■ Vti j = -Vp - /9V$ - pV$e^t + /. 



(1) 



(2) 



Here p is the density, u is the velocity field, $ is the gravi- 
tational potential due to the disc, ^ext is the potential due 
to external bodies, p is the pressure and / is the viscous 
force. The gravitational potential due to the disc satisfies 
Poisson's equation 



(3) 
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which yields the integral expression 

p{r\ip\z'ydr'dip'dz' 



-G 
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(4) 



where V is the domain where p is non zero. This coincides 
with the disc domain when this is not separated from exter- 
nal material. 

For analytic discussion we consider fluids with a gen- 
eral barotropic equation of state for which p — p(p) and 



dp /dp 



with Cs being the local sound speed; for nu- 



merical calculations a locally isothermal equation of state is 
adopted, p = c^p, where Cs{r,z) is a specified function of r 
and z. 



2.1 Razor thin limit 

To obtain the limiting case of a razor thin disc, we assume no 
interior vertical motions and that the radial and azimuthal 
velocity components are independent of z. We then inte- 
grate over z so that p is replaced by the surface density S in 
equation ([T]) and equation ([2]) while p becomes a vertically 
integrated pressure, which in the barotropic case is assumed 
to be a function only of S with dp/dYl — c^. The z depen- 
dence of the potentials is neglected so that we now have the 
governing equations 

^+V-(uE)=0, (5) 

and 



du 
~dt 



+ u • Vu 



-Vp - - SV^ext + /, 



(6) 



where u = {ur^u^) is now a two dimensional velocity field 
and / is now the two dimensional viscous force, charac- 
teris ed by a unif orm kinematic viscosity v in our models 
(see lMassetll2002l , for details). The midplane disc potential 
is given by 



y^r^ + r''^ — 2rr' cos — i^') + 



rdrd(p\ (7) 



where we have introduced a softening length Cg — egoH(r'), 
with Cgo being a dimensionless constant and a putative semi- 
thickness for razor thin discs defined through H(r) = hr — 
Cs/^k where Qk = \/ GM^ / is the Keplerian rotation rate. 
The dimensionless constant disc aspect ratio is h. 

The gravitational potential due to external bodies com- 
prising the central star and an embedded planet is 

GM, GMr, 



+ r 

Id 

^ GM, 



Y^r^ + r| - 2rrp cos {ip - (pp) + 
— cos [cp — (f )r dr dcp 

-r cos (cp — ifp). 



(8) 



Here Mp is the fixed mass of the embedded planet with 
cylindrical coordinates {rp{t), (pp{t)). The associated grav- 
itational potential is softened with softening length Cp = 
epoH{rp), where Cpo is a dimensionless constant. The last 
two terms in the expression for ^ext which give the indirect 
potential account for the forces due to the disc and embed- 
ded planet acting on the central star. They occur because 
we adopt a non-inertial frame of reference. 



2.2 Model setup 

We describe specific disc and planet models used in numeri- 
cal calculations, which also motivates the analytical discus- 
sion below. We adopt units such that G = M* = 1 and the 
inner boundary radius of the disc, r^, is unity. The Keplerian 
orbital period at this radius r = = 1 is then P(l) = 27r. 
The disc occupies r = [n,ro] = [1, 10]. The disc is taken to 
have a locally isothermal equation of state, p = CgS, where 
Cs = rhQk- We remark that softening is introduced to take 
account of the vertical thickness of the disc. Typically, soft- 
enin ^ length to s e mi-thickness ratios of 0.3 — 0.6 are used 
fe.g. lMassetll2002l : iBaruteau Massetll2008h . Thus fiducial 
values for the constants Cpo, ^go and /i, of 0.6, 0.3 and 0.05 
are respectively adopted. The initial surface density profile 
is taken to be given by 



-3/2 



1 



(9) 



(see'Armitage & Hansen'l999") where Hi = H{ri). This form 
for S, which vanishes smoothly a distance Hi inside the disc 
inner boundary, has been introduced to ensure that both the 
gravitational force and the pressure gradient in hydrostatic 
equilibrium at the disc inner boundary give minor contribu- 
tions and vary smoothly there. The constant surface density 
scale So is adjusted so that (5(ro), where 



(10) 



ttGS 7rr2S(r) 

is the Toomre Q parameter, evaluated in the limit of a thin 
Keplerian disc for which the epicycle frequency n — Q^k^ 
takes on a specified value Qo. Thus Qo parametrises the 
models. 

The initial azimuthal velocity is found by assuming 
hydrostatic equilibrium in which the centrifugal force bal- 
ances forces due to stellar gravity and the disc's self-gravity 
and pressure gradient. Thus the initial azimuthal velocity is 
given by 



2 r dp GM^ d^ 
1^ dr r dr 



(11) 



For the local isothermal equation of state we have adopted, 
the contribution due to the pressure is given by 



T^dp_ _ 2 



5 

-2 + 



ry/n{r + Hi) 



-3/2 



2 [l - Vn/(r + Hi)\ J ■ 



(12) 



At r = for h — 0.05, this is 2OC5 which is approx- 
imately 5% of the square of the local Keplerian speed so 
that the contribution of the pressure force is indeed minor 
when compared to that arising from the gravity of the cen- 
tral star. 

The initial radial velocity is set to Ur — 3z^/r, corre- 
sponding to the initial radial velocity velocity of an one- 
dimensional Keplerian accretion disc with S (x r~^^'^ and 
uniform kinematic viscosity. The disc is evolved for a time 
^ 280P(ri) before introducing the planet. 

When a planet of mass Mp is introduced, it is inserted 
in a circular orbit, under the gravity of both the central star 
and disc, at a distance rp — rp{t — 0) from the central star. 
We quote time in units of the Keplerian orbital period at the 
planet's initial radius, which is given by Po = 2ti /Q.k{Tp{t — 
0)). If the planet is held on fixed circular orbit then rp(t) — 
rp{t = 0). 
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3 NUMERICAL SIMULATIONS 

In order to motivate and facilitate our analytic discussion, 
linear analysis and interpretation of the instabilities we 
find in a self- gravitating disc with a surface density gap or 
dip induced by a massive planet, we here provide a brief 
demonstration of their existence. We show global instabili- 
ties associated with a surface density depression made self- 
consistently by a giant planet in a fixed circular orbit by 
means of numerical simulations. Details of the numerical 
approach are given in ^ where additional hydrodynamic 
simulations of this type as well as others, where the planet 
is allowed to migrate, will be investigated more fully. 



3.1 The nature of edge modes at the edges of disc 
surface density depressions induced by 
interaction with giant planets 

Here, we describe results obtained for a disc with the ini- 
tial density profile specified above in which there is an em- 
bedded planet with mass Mp = 3 X 10"^M*. This corre- 
sponds to a Saturn mass planet when M* = IMq. The 
disc model is Qo = 1-5, corresponding to total disc mass of 
Md — 0.063M*. The uniform kinematic viscosity was taken 
to be = 10~^ in code units. The planet was introduced 
at radius = 5 at t = 25P(rp) = 25Po. Its gravitational 
potential was then ramped up over a time interval lOPo- For 
this simulation it was held on a fixed circular orbit. Without 
a perturber we expect, and find, the disc to be gravitation- 
ally stable because we have Q{rp) = 2.62 and near the outer 
boundary Q ^ 1.5. 

The profile of the gap opened by the planet is shown 
in Fig. [T] The azimuthally averaged surface density (S)c^, 
vortensity {r])^ where 77 = Ai:^/2QS, and Toomre Q value 
are plotted. The latter is calculated using the azimuthally 
averaged epicycle frequency. We remark that for an axisym- 
metric disc, the Toomre parameter is proportional to the 
product of the ratio of the rotation frequency Q to the epicy- 
cle frequency k and the vortensity 77. It is seen that there 
is a surface density depression, referred to a gap, associated 
with a decrease of surface density by ^ 20% relative to the 
unperturbed disc. 

It has been found that extrema in the vorten sity are as- 
sociated with instability (Papaloizou & Lin 1989). For typi- 
cal disc models structured by disc-planet interactions, such 
as those illustrated in Fig.[T] local maxima/minima in Q and 
r] have been found to approximately coincide. The neigh- 
bourhoods of the inner and outer gap edges both contain 
maxima and minima of Q and 77. 

In the absence of a planet, Q smoothly decreases out- 
wards. In the case when a planet is present, disc-planet in- 
teraction results in signific ant vort ensity generation as ma- 
terial flows through shocks (|Lin Papaloizou 2010), leading 
to vortensity maxima. The resulting Q and (rj)^ profiles are 
very similar (since Q = (cs/7rG')(2Q?7/S)"^/^). Fig. [T] shows 
that the planet-modified Q profile may exhibit a range of 
behaviours close to r^. Vortensity diffusion, e.g., due to sig- 
nificantly larger viscosities than those considered here, could 
render the Q profile to be uniform. In such cases, unsta- 
ble modes associated with vortensity maxima cannot be set 
up. The unperturbed Toomre-Q profile (prior to the planet 
introduction) may also play a role. The Toomre-Q profile 
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Figure 1. Gap profile produced by a Saturn-mass planet embed- 
ded in a disc with Qo = 1.5. The azimuthally averaged surface 
density (dotted line), vortensity (dashed line) and Toomre Q pa- 
rameter (solid line) are plotted. Note that the Q profile shows 
a maximum and minimum associated with each gap edge. The 
planet is in a fixed circular orbit located at r = 5. 



perturbed by the planet is reminiscent of its unperturbed 
value, without planet. If the latter happens to be approxi- 
mately uniform, the net impact of edge modes on the planet 
torque could be negligible. In our models, the background 
Toomre Q decreases with radius, so we expect that the outer 
gap edge is more gravitationally unstable than the inner gap 
edge. 

Returning to the present case. Fig. [T] shows that in the 
neighbourhood of the outer gap edge, the maximum and 
minimum values of Q occur at r = 5.45 and r = 5.75. and are 
Q — 4.2 and Q — 1.75 respectively. In the region of the inner 
gap edge, the maximum and minimum values of Q occur at 
r = 4.55 and r = 4.25. These are Q = 4.75 and Q = 2.55 
respectively. The extrema are separated by lAH and I.IH 
in the inner and outer gap edge regions respectively. Thus 
the characteristic width of the gap edge is the local scale- 
height. On average Q ^ 2 for r > 6. Since Q > 1 everywhere, 
the disc is stable against local axisymmetric perturbations. 

The simulation shows that both gap edges become un- 
stable. The surface density contours at t = 50Po are shown 
in Fig. [21 We identify a mode with m = 3 and relative density 
perturbation AS/S ~ 0.4 associated with the inner edge 
and a mode with m — 2 and AS/S :^ 0.9 associated with 
the outer edge. The modes have become nonlinear, with the 
development of shocks, on dynamical time-scales. Spirals do 
not extend across the gap indicating effective gap opening 
by the planet and disc self-gravity is not strong enough to 
connect the spiral modes on either side of Vp. 

Plots of the radial dependence of the pattern rotation 
speed, r^pat, obtained from the m = 2 component of the 
Fourier transform of the surface density, together with the 
azimuthally averaged values of Q and Q.^hi/2 are presented 
in Fig. [21 We focus on the mode associated with the outer 
gap edge because the m — 2 spiral mode has more than twice 
the relative surface density amplitude compared to the inner 
spiral mode. 

By necessity, the Fourier transform includes features 
due to the planetary wakes and inner disc spiral modes 
as well as the dominant outer disc mode. On average. 
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^pat 0.08 El- There is a co-rotation point in the outer 
disc at r = 5.5 where Q = Qpat 0.078, i.e. at the lo- 
cal vortensity (surface density) maximum (minimum) of the 
basic state, which is just within the gap (see Fig. [T]). This is 
consistent with pattern speeds obtained from linear calcu- 
lations presented later. Fig. [2] indicates another co-rotation 
point at r ^ 4.7 but this is due to inclusion of features as- 
sociated with the planetary wakes. 

We refer to unstable spiral modes identified here as 
edge modes. Their properties c an be compared to thos e of 
groove mod es in particle discs ([Sellwood Kahn"l 99lh or 
fluid discs (jMeschiari k LaughlinI l2008l ) . ISellwood KahnI 
describe groove modes as gravitational instabilities associ- 
ated with the coupling of disturbances associated with two 
edges across the gap between them. However, our analytic 
description of the edge mode instability is a coupling be- 
tween a disturbance at a single gap edge and the distur- 
bance it excites in the adjoining smooth disc away from co- 
rotation. The other gap edge plays no role. Although both 
types of mode are associated with vortensity maxima, the 
groove modes described above would have co-rotation at the 
gap centre midway between the edges, whereas in our case 
co-rotation is at the gap edge. 

This is significant because in our model there will be dif- 
ferential rotation between the spiral pattern and the planet. 
In addition, unstable modes on the inner and outer gap edges 
need not have the same m (Fig. [2] shows m = 3 on the inner 
edge and m = 2 on the outer edge). If there was a groove 
mode with co-rotation at the gap centre, then the spiral pat- 
tern would co-rotate with the planet and the coupled edges 
must have disturbances dominated by the same value of m. 



3.2 Comparison to vortex instabilities 

We have demonstrated a global instability in self- gravitating 
disc-planet systems with a gap. It is interesting to com- 
pare these spiral modes to the well-known localised vortex- 
forming instabilities that occur in non- self- gravitati ng or 
weakly self-sravitating discs near gap edges (|Li et al .11 20091 : 
iLin k PaDaloizoull201Ql ). However, the vortex instability re- 
quires low viscosity. For comparison purposes, here we use a 
kinematic viscosity of = 10~^ (or an a viscosity parameter 
0(10-^)). This value of u is an order of magnitude smaller 
than what has been typically adopted for protoplanetary 
disc models. We compare the behaviour of disc models with 
Qo = 1.5 and Qo = 4.0, the former having strong self-gravity 
and the latter weak self-gravity. 

Fig. [3] shows two types of instabilities depending on disc 
mass. The lower mass disc with Qo = 4 develops six vortices 
localised in the vicinity of the outer gap edge. The more 
massive disc with = 1.5 develops edge modes of the type 
identified in the earlier run with v = 10~^. Here, the m = 
3 spiral mode is favoured because of the smaller viscosity 
coeflficient used. 

Vortex modes are said to be local because instabil- 
ity is ass ociated with fiow in the vicinity of co-rotation 
([Lovelace et al.ll 19991 ) and does not require the excitation of 



We have explicitly checked that this is the pattern speed by 
measuring the angle through which the spiral pattern rotates in 
a given time interval. 
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Figure 2. Relative surface density perturbation for Qo = 1.5 at 
t = 50Po (top) and the radial form of the pattern speed found 
from the m = 2 component of the Fourier transform of the surface 
density (bottom, solid). Note that this mode is associated with 
the outer disc and, apart from within the gap, this mode appears 
to have a stable pattern speed of Qpat ~ 0.078 corresponding to 
a corotation point at r = 5.5. Plots of Q and Q ^ hi/2 are also 
given, (bottom, dashed). 



waves. Edge modes are said to be global because instability 
requires interaction between an edge disturbance and waves 
launched at Lindblad resonances, away from co-rotation, in 
the smooth part of t he disc. Vortices perturb th e disc even 
without self-gravity ([Paarde kooper et al.l I2OIOI ) . Although 
the vortex mode in Fig. [3] shows significant wave pertur- 
bations in the outer parts of the disc, these should be seen 
as a consequence of unstable vortices forming around co- 
rotation, rather than the cause of a linear instability. 

We found the vortex modes have co-rotation close to 
local vortensity minima in the undisturbed gap profile. The 
vortex modes are localised with high m being dominant 
(m > 5). On the other hand the spiral modes found here 
are global with low m < 5. Edge modes make the gap less 
identifiable than vortex modes do because the former, with 
corotation closer to r^, protrudes the gap edge more. While 
the edge mode only takes a few Pq to become non-linear 
with spiral shocks attaining comparable amplitudes to the 
planetary wakes, the vortex mode takes a significantly longer 
time (the plot for Qo = 4 is 30Po after gap formation). 
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Figure 3. Relative surface density perturbation for the Qo = 1.5 
disc (left) and the Qo = 4.0 disc (right) with u = IQ-^ . These 
plots show that the gap opened by a Saturn-mass planet supports 
vortex instabilities in low mass discs and spiral instabilities in 
sufficiently massive discs. 



4 ANALYTICAL DISCUSSION OF EDGE 
MODES 

The fiducial disc simulated above is massive with Md ^ 
0.06M* but it is still stable against gravitational instabilities 
in the sense that Q ^ 1.5 everywhere and the disc without 
an embedded planet, which has a smoothly varying surface 
density profile, does not exhibit the spontaneous growth of 
spiral instabilities. When a planet of Saturn's mass is in- 
serted (Mp = 3 X 10~^M*), it is expected to only open a 
partial gap (^ 30% deficit in surface density), but even this 
is sufficient to trigger a m = 2 spiral instability. 

The fiducial case above thus suggests spiral modes may 
develop under conditions that are no t as e xtreme as those 
considered by iMeschiari LaughhnI (|2008l ). They used a 
prescribed gap profile with a gap depth of 90% relative 
to the background, which is three times deeper than in 
our models. Their gap corresponds to Mp = 0.002M*, al- 
though a planet potential was not explicitly included. In 
both their model and ours, the gap width is ^ 2rh where 
rh = {Mp/3M^y^^rp is the Hill radius, but since they effec- 
tively used a two Jupiter mass planet, at the same rp their 
gap is about 1.9 times wider than ours. 

We devote this section and the next to a theoretical dis- 
cussion of spiral modes associated with planetary gap edges. 
This work is based on the an analysis of the governing equa- 
tions for linear perturbations. As angular momentum bal- 
ance is important for enabling small perturbations to grow 
unstably, we begin by formulating the conservation of angu- 
lar momentum for linear perturbations. 



4.1 The conservation of angular momentum for a 
perturbed disc 

We derive a conservation law for the angular momentum as- 
sociated with the perturbations of a disc that enables the 
angular momentum density and flux to be identified within 
the framework of linear perturbation theory. The behaviour 
of these quantities is found to be important for indicating 
the nature of the angular momentum balance in a system 
with a neutral or weakly growing normal mode and how pos- 
itive (negative) angular momentum fluxes associated with 
wave losses may drive the instability of a disturbance that 
decreases (increases) the local angular momentum density. 



4.2 Barotropic discs 

We begin by writing down the linearised equation of mo- 
tion in three dimensions for the Lagrangian displacement 
^ = (^rj^ifj^z) for a differentially rotating fluid with a 
barotropic equation of state and with self-gravity included 
(see eg. lLvnden- Bell Ostrikei]ll967l : lLin et al.ll 19931 ) in the 
form 



Dt^ 
where 



Dt 



5^ 



p 



(13) 



(14) 



Here we adopt a cylindrical polar coordinate system (r, cp, z), 
perturbations to quantities are denoted with a prime while 
unprimed quantities refer to the background state, the op- 
erator D/Dt = d/dt -\- Qd/dcp is the convective derivative 
following the unperturbed motion and k is the unit vector 
in the z direction which is normal to the disc mid-plane. For 
perturbations that depend on cp through a factor exp(im(/?), 
m being the azimuthal mode number, assumed positive, that 
we consider, the operator D/Dt reduces to the operator 
{d/dt-\-imQ) . In addition we recall that for a barotropic 
equation of state, Q = Q(r) depends only on the radial co- 
ordinate. Then ([13)) becomes 

^ + 2Qk A ^ + 2imQ^ + 2imQ^fc A | 
ot^ at at 



The perturbation to the density is given by 

p' = - V • {pO- 



(15) 



(16) 



The perturbation to the disc's gravitational potential is 
given by linearising Poisson's equation to give 

V^$' = AttGp. (17) 

We have also added an external potential perturbation $ext 
to assist with the identification of angular momentum flux 
later on. 

By taking the scalar product of ([15]) with multi- 
plying by the background density, p, and taking the imag- 
inary part, after making use of ([T6)) and ([T7)) we obtain a 
conservation law that expresses the conservation of angular 
momentum for the perturbations in the form. 



%- + V • (Fa + Fg + Fe> 
at 

where 



pj^- ^Im . g + Qk . ^ A r + imm' 



2 



(18) 

(19) 
(20) 
(21) 
(22) 



Fext 

and 

T = ^Im (<f Ltp'*) . (23) 
Here pj is the angular momentum density and the angular 
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momentum flux is split into three contributions, Fa being 
proportional to the Lagrangian displacement, is the advec- 
tive angular momentum flux, Fg is the flux associated with 
the perturbed gravitational stresses and Fext is the flux as- 
sociated with the external potential which is inactive for free 
perturbations. 

The quantity T is a torque density associated with the 
external potential as can be seen from that fact that the real 
torque integrated over azimuth and divided by 27v is 

-ir (^) '^=^= > ^"'-'''^^ • 

(24) 

This justifies the scalings used for the angular momentum 
density and fluxes. 



An expression of the form ([T8|) may be derived but now 
the torque density T takes the form 

T = ^ [Im - Im (p'C ■ V(cJ))] . (29) 

When Cs is constant corresponding to a strictly isothermal 
equation of state, the fluid is again barotropic and previous 
expressions recovered. Equation [29l contains terms in addi- 
tion to external contributions. In general these imply the 
possibility of an exchange of angular momentum between 
perturbations and the background. A very similar discus- 
sion applies to the razor thin limit. We now go on to apply 
the above analysis to determining the angular momentum 
balance for the linear perturbations of razor thin discs in 
more detail and thus determine properties of and conditions 
for unstable normal modes. 



4.3 The conservation of angular momentum in 
the two dimensional razor thin disc limit 

The form of the conservation law for a razor thin disc is 
obtained by integrating equation ([T8|) over the vertical coor- 
dinate. For terms oc p the integrand is non zero only within 
the disc. ^ and the gravitational potentials may be assumed 
to depend only on r and cp. The effect of the integration is 
simply to replace p and p' by the surface density S and its 
perturbation respectively. The fluxes become vertically 
integrated fluxes. For the term associated with the gravita- 
tional stresses, assuming the potential perturbation vanishes 
at large distances, we have 

IdirFcr) 



L 



dr 



where 



(25) 



(26) 



Thus for a two dimensional razor thin disc the conserva- 
tion law is of the same form as (|18p but with pj. Fa, and 
¥ext given by {[19]), (|2Q)) and ([22]) with p replaced by S re- 
spectively. As the vertical direction has been integrated over 
only the radial components of the fluxes contribute. From 
the above analysis the vertically integrated angular momen- 
tum flux due to gravitational stresses becomes 



ForV, 



(27) 



with Fcr given by ()26p . Accordingly this term, unlike the 
others involves an integration over the vertical direction. We 
shall obtain explicit expressions for the fluxes in terms of the 
perturbation S' — T^' c?s /^-\-^' appropriate to the razor thin 
disc below. 



4.4 The conservation of angular momentum for a 
disc with a locally isothermal equation of state 

It is possible to repeat the analysis of section 14.11 when the 
disc has a locally isothermal equation of state. In this case 
P — pc?s^ where Cs is a prescribed function of position and 
the linearised equation of motion ([13]) is modified to read 



-VS' + ^V(c^)-V<,t. 
P 

(28) 



4.5 Properties of the linear modes of a razor thin 
disc 

We assume linear perturbations for which the Lp and t de- 
pendence is through a factor of the form g^C^^+^v') ^here 
a is the complex eigenfrequency and m is the azimuthal 
mode-number. From now on this factor is taken as read. 
Linearising the hydrodynamic equations ([5]) and ([6]) for the 
inviscid case and ignoring external potentials, we obtain 



1 d{T.ru'r) imllu'^ 



r dr 
dS' ^ T.' del 
dr S dr 
imS' 



(30) 
(31) 
(32) 



where a = cr + mQ is the shifted mode frequency. Using equa- 
tions ()3ip and ()32l) to eliminate u'^ and u' in e quation (|3Qp 
we obtain (see eg. IPapaloizou k, Savoniielll99lh 



dr 



r^ 
D 

2mQai; 



dS' 2mVLaS' 



dS' 2mQaS' 



rS^ dcV 
D dr 



mS' 



dr 



+ ■ 



cr dr 
2mQ.T.' del 
Dra dr 



(33) 



where D = — a^, and recall 77 = k'^ /2QTi is the vortensity. 
This form shows that a corotation singularity where a — 
can be avoided if corotation corresponds to a vortensity 
extremum. 

The above expressions may be applied to a disc, ei- 
ther with a locally isothermal equation of state, or with a 
barotropic equation of state, where the integrated pressure 
P = P(S) and the soundspeed cl = dP/dT^. To obtain the 
latter case , the quantity dc?s/dr is simply replaced by zero. 
We remark that inclusion of additional terms dependent on 
this has been found numerically to produce only a slight 
modification of the discussion that applies when they are 
neglected. This is because varies slowly compared to ei- 
ther the linear perturbations or the surface density in the 
vicinity of the edge, thus we shall neglect (r / c?s)d(?s / dr from 
now on. We also recall that S' — T^' cl /T^+^' and the gravita- 
tional potential is given by the vertically integrated Poisson 



8 Lin and Papaloizou 



integral as 

<t>' = s' 



Km{r,r')T^'{r')r'dr' , 



where 



f 

Jo 



cos{m(p)d(p 



(34) 



(35) 



■ 2rr' cos (cp) + 

Here the domain of integration T> is that of the disc pro- 
vided there is no external material. However, in a situation 
where density waves propagate beyond the disc boundaries 
T> either should be formally extended or the form of Km 
changed to properly reflect the boundary conditions. We re- 
mark that ()34p enables to be determined in terms of S' 
and if this is use d in ([33l) a single eigenvalue equation for 
results (eg. Papa loizou Savonije 1991). 

We shall now consider the angular momentum flux bal- 
ance for normal modes. For simplicity we shall consider 
the barotropic case. The analytic discussion concerns a disc 
model with one boundary being a sharp edge with arbi- 
trary propagation of density waves directed away from this 
boundary being allowed. When applied to a gap opened by 
a planet in a global disc, the analytic model corresponds 
the section of the disc from the inner disc boundary to the 
inner gap edge or the section from the outer gap edge to 
the outer disc boundary. These two sections are assumed to 
be decoupled in the analytic discussion of stability but not 
in the numerical one. In addition, although it produces the 
background profile, the planet is assumed to have no effect 
on the linear stability analysis discussed here, but it is fully 
incorporated in nonlinear simulations. The correspondence 
between the various approaches adopted indicates the above 
approximations are reasonable. The discussion of the outer 
and inner disc sections is very similar, accordingly these are 
considered below together. 



4.6 Angular momentum flux balance for normal 
modes 

The vertically integrated angular momentum flux associated 
with perturbations can be related to the background vorten- 
sity profile. For this discussion we use the governing equation 
in th e form given by eg nation (|33|) ("Papaloizou & Savonije 
Il99lh . We assume a barotropic disc model and neglect terms 
involving dcs/dr. Multiplying this equation by S'* and inte- 
grating, we get 



HS', 



nr 
Jri 

Jr, D \dr 

i: 



rT.'S"'dr± 



dS' 2mQd- 



+ 



dS' 2mVLa 



+ ■ 



dr rn- 

202^-11 |2 

dr + 



S' 



nr 



dr 

dS^ , 2mQa 
dr 

m|5f 



s' s"- 



+ ■ 



S"" dr 



dr 



-\dr ■ 



0, (36) 



where, here and below, the positive sign alternative appUes 
to the case where the disc domain has an inner sharp edge 
where r — ri and then rb = To, the outer boundary radius. 
The negative sign alternative applies to the case where the 
disc domain has an outer sharp edge where r — ro and then 
rh — Ti^ the inner boundary radius. In both cases waves may 
propagate through r = rb. Note also that there is no contri- 
bution from edge boundary terms as the surface density is 
assumed to be negligible there. 



We may now express the terms in the above integral 
relation, which for convenience we have taken to define a 
functional of S' that depends on a, in terms of quantities 
related to the transport of angular momentum by taking its 
imaginary part. We shall begin by assuming that a is real 
or, more precisely, that we are in the limit that marginal 
stability is approached. 

By making use of equation (|25]) we find that 

2 



Im 



■ rro 

.J ri 



rT.'S^'dr 



= Im 



/ro 
■ ■ i 



rS' 



dr 



Im 



nr 



rT,'<t>'*dr 



(37) 



We recall that Fcr is the vertically integrated angular mo- 
mentum fiux transported by gravitational stresses. Note that 
Fcr at the sharp edge is ignored. Equation l25[ together with 
the requirement that Fcr is regular for r ^ and vanishes 
more rapidly than 1/r for r ^ oo, implies that Fcr is zero at 
the sharp edge separating the disc domain and the exterior 
(assumed) vacuum or very low density region. 

In addition we remark that from equations pTI ) and 
(|32)) . the radial Lagrangian displacement is given by 



c — r _ 

Sr — — — 
ICT 



1 

D 



dS' ^ 2mS'VL 
dr ra 



(38) 



From this it follows that the imaginary part of the sec- 
ond term on the right hand side of equation (|36)) is 
— it Im [rTi^rS'*]\ri, = —^{2/m)[rFAr]\ri,- Using the above, 
taking the imaginary part of ()36|) yields the remarkably sim- 
ple expression 



± [r{FGr + FAr 



- m\S^d_ 
7 dr 



-]dr 
V 



(39) 

Note that we have been assuming that a is real and that 
the mode is at marginal stability. Then the right hand side of 
(|39|) is apparently real. However, the integrand is potentially 
singular at corotation where a — 0. Thus the approach to 
marginal stability has to be taken with care. 

Setting a = an — ij, where an, being the real part of a, 
defines the corotation radius, rc, through Q{rc) = —aR/m 
and the growth rate as 7, with —7 being the imaginary part 
of a. Marginal stability can be approached by assuming that 
7 has a vanishingly small magnitude but is positive. Then 
we can use the Landau prescription and substitute l/a by 
V{l/d-) +i7r^(a), where V indicates that the principal value 
of the integral is to be taken and 6 denotes Dirac's delta 
function. Adopting this (|39p becomes 

\sr c 



[riFcr + FAr 



±- 



\dn/dr\ dr 



(40) 



The above relation can be viewed as stating that at marginal 
stability, either corotation is at a vortensity extremum and 
both terms in equation (|4Qp vanish, or angular momen- 
tum losses as a result of waves passing through r = rb, 
are balanced by torques exerted at corotation. The latter 
torqu e is proportion al to the gradient of the vortensity , 77 
(Goldreich Sz Tremaine 1979) and recall that r]~^ scales with 
E. 

The propagating waves quite generally carry negative 
angular momentum when located in an inner disc and pos- 
itive angular momentum when located in an outer disc 
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(jGoldreich Tremainelll979l l. Thus a balance may be pos- 
sible between angular momentum losses resulting from the 
propagation of these waves out of the system and angular 
momentum gained or lost by an edge disturbance as ex- 
pressed by equation (|4Q|) . For an edge disturbance in a region 
of increasing (decreasing) surface density, such as an inner 
(outer) edge to an outer (inner) disc, the edge disturbance 
loses (gains) angular momentum. 



In practice, an inner (outer) edge disturbance may be 
associated with a positive (negative) surface density slope 
near a vortensity maximum or near a vortensity minimum 
occurring as a result of the variation of both and S. The 
former case is associated with discs for which self-gravity 
is important and spiral density waves are readily excited. 
The latter is associated with weakly or non self- gravitating 
discs and is associated with vortex formation at gap edges 
as has already been dis cussed by several authors (see eg. 
iLin Papaloizoul (|2010l l and references therein). 



So far we have not discussed the boundary condition 
at r = (the non-sharp boundary). One possibility is that 
this corresponds to stipulating outgoing waves or a radiation 
condition, so that the fluxes in (|4Qp are non zero. Another 
possibility is to have a surface density taper to zero which 
would mean wave reflection and removal of the boundary 
flux terms. In fact we expect that the issue of stability is 
not sensitive to this. We remark that in the case of non self- 
gravitating discs and the low m modes of interest, the re- 
gions away from the edge are evanescent and amplitudes are 
small there, the disturbance being localised in the vicinity 
of the edge. For self-gravitating discs, as we indicate below, 
instability can appear because of an unstable interaction 
between outwardly propagating positive {inwardly propagat- 
ing negative) angular momentum density waves and a nega- 
tive (positive) angular momentum edge disturbance localised 
near an inner (outer) disc gap edge. As long as these waves 
are excited it should not matter whether they are reflected 
or transmitted at large distance, as long as their angular 
momentum density remains in the wider disc and is not fed 
back to the exciting edge disturbance. 



For these reasons and also simplicity we shall assume 
that at the boundary where r = rb, either the surface den- 
sity has tapered to zero, or there is reflection such that the 
boundary terms in expressions like (|36p may be dropped. 
Equation ()40|) then simply states that marginal stability oc- 
curs when corotation is at a vortensity extremum. We find 
that the vortensity maximum case is associated with an 
instability in self- gravitating discs that is associated with 
strong spiral waves. On the other hand the vortensity min- 
imum case is associated with the vortex- forming instability 
in non-s elf gravitating disc s that has been previously stud- 
ied (see iLin Papaloizoi] (|201Q| ) and references therein) 
and will not be d i scusse d further in this paper ( but see 
iLin Papaloi^ dioiTI l for a discussion of the effect of 
weak self-gravity on vortex modes). 



5 DESCRIPTION OF THE EDGE 

DISTURBANCE ASSOCIATED WITH A 
VORTENSITY MAXIMUM IN A 
SELF-GRAVITATING DISC 

We now focus on the description of the instability in a self- 
gravitating disc as being due to a disturbance associated 
with the edge causing the excitation of spiral waves that 
propagate away from it resulting in destabilisation. This oc- 
curs because the emitted waves carry away angular momen- 
tum which has the opposite sign to that associated with the 
edge disturbance ( see section and equation (|40)) ). Ac- 
cordingly the excitation process is expected to lead to the 
growth of this disturbance. 

We thus consider the disturbance to be localised in the 
vicinity of the edge and the potential perturbation it pro- 
duces to excite spiral density waves in the bulk of the disc. 
Our numerical calculations show that edge dominated modes 
of this type occur for low m and are dominant in the non- 
linear regime. We assume the mode is weakly growing cor- 
responding to the back reaction of the waves on the edge 
disturbance being weak. Thus in the first instance we calcu- 
late a neutral edge disturbance and then calculate the wave 
emission as a perturbation. We emphasise again that an im- 
portant feature is that corotation is at a vortensity maximum 
(in contrast to the non self- gravitating case for which coro- 
tation is located at a vortensity minimum). We show that 
such modes require self-gravity. In an average sense they re- 
quire a sufficiently small Q value and so cannot occur in the 
non self-gravitating limit (Q oo). 



5.1 Neutral edge disturbances with corotation at 
a vortensity maximum 

We start from equation ()33p for S^ We assume that only 
the third term on the right hand side need be retained. This 
is because this term should dominate near corotation in the 
presence of large vortensity gradients that are presumed to 
occur near the edge and the disturbance is assumed localised 
there (see also the discussion of groove modes in collision- 
less particle discs of lSellwood KahnI (|l99lh , where similar 
assumptions are made). Thus we have 



ruj dr 



(41) 



where cj — a/m. The associated gravitational potential is 
given by 



-G 



Km(r,r') 



S'ir') 



dr' . 



(42) 



From the relation S ■ 
gral equation 



u(r') dr' \_ri(r') 
S^Cs/S + ^' we thus obtain the inte- 



S'{r) 



rHuj dr 



z(r,r)-^ 



uo(r') dr' 



r](r') 



dr' . 



(43) 



For a disturbance dominated by self-gravity, pressure should 
be negligible. Under such an approximation we have S' ^ ^' . 
Equation ()41|) then implies that to obtain a negative (posi- 
tive) potential perturbation for a positive (negative) surface 
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density perturbation, we need sgn(cj) = sgn(dri/dr). At co- 
rotation where uj — and there is a vortensity extremum, 
this requirement becomes sgn{dCu/dr) = sgn{dQ/dr) = 
sgn(d^r]/dr^). Since typical rotation profiles have dQ/dr = 
< 0, the physical requirement that potential and sur- 
face density perturbations have opposite signs implies that 
d^rj/dr^ < at co-rotation, i.e. the vortensity is a maximum. 
Thus our discussion applies to vortensity maxima only. 
For the class of vortensity profiles for which sgn(dri/dr) — 
sgn(c<;), we demonstrate below that the integral equation 
(|43|) may be transformed into a Fredholm integral equation 
with symmetric kernel. 

Noting that corotation is located at a vortensity maxi- 
mum, the requirement above imply the vortensity increases 
(decreases) interior (exterior) to corotation ( we comment 
that the discussion may also be extended to the case when 
that is true with d^r]/dr^ vanishing at corotation). As we 
expect the disturbance to be localised around corotation, it 
is reasonable to assume this holds and that we may contract 
the integration domain (ri,ro) to exclude regions which do 
not conform without significantly affecting the problem. 

Proceeding in this way we introduce the function H 
such that 



where 
Z{r) 

Defining the new symmetric kernel 7l{r,r') as 
7^(r, r') =GK,^ir,r')yir)yir') 





[-)-] 


-1/2 






dr 






rEcj dr 


m 



-1/2 



where 



y(r) 



dr \r] J uj 



11/2 



1 



rScj dr \r] 



-1/2 



we obtain the integral equation 

pro 

n= n[r,r')H[r')dr'. 

In other words, 1-i is required to be the solution to 
Xn= n{r,r')H[r')dr' , 



(44) 



(45) 



(46) 



(47) 



(48) 



(49) 



with A = 1. However, Eq.[49]is just a standard homogeneous 
Fredholm integral equation of the second kind. Thus the 
existence of a nontrivial solution (1-L ^ 0) implies 1-L is an 
eigenfunction of the kernel IZ with unit eigenvalue. 

The kernel IZ is positive and symmetric. Accordingly 
the Fredholm equation above has the property that the max- 
imum eigenvalue A > and is the maximum of the quantity 
A given by 



A = 



/;;/;; 7^(r,r07^(/)7^*(r)drd/ 

!:^mr)\^dr 



(50) 



over all H,. Thus, if it is possible to show that we must always 
have max(A) < 1, then it is not possible to have a non-trivial 
solution corresponding to a neutral mode. Upon application 
of the Cauchy-Schwarz inequality (twice), one can deduce 



/To 



\lZ(r^r')^ drdr' — A^^st ^ max(A^). 



(51) 



Thus the existence of a neutral mode of this type is not 
possible if A^st is less than unity. 

To relate the necessary condition for mode existence to 
physical quantities, we can estimate h^st- The Poisson kernel 
Km{r, r') is largest at r = r^ Other factors in the numerator 
of IZ involve the factor 1/cj. We assume these factors have 
their largest contribution at corotation where r — rc and 
uj — 0. Hence, A^st Ir° Sr° \7l(rc^rc)\^ drdr' . Assuming 
the edge region has width of order Lc, taken to be much less 
than the local radius but not less than the local scale- height, 
we estimate 



Aest '■^\R(rc,rc)\Lc 

- GKm{rc,rc) 



(52) 



All quantities are evaluated at co-rotation and the double 
prime denotes d^ /dr^ . Because the edge is thin by assump- 
tion, we can further approximate the Poisson kernel by 



Km{r,r') 



'rr' 



-.Ko 



/rr' 



(53) 



where Kq is the modified Bessel function of the second kind 
of order zero and recall that eg is the softening length. If 



pressure effects are negligible 
. 2GKo{meg/rc) 

Aest ~ 



{cs being small), we obtain 



1 









(54) 



Since 77 = hi^ /2QT^, the requirement that Aest exceeds unity, 
leads us to expect that modes of the type we have been 
discussing can exist only for sufficiently large surface density 
scales. We have approximately 



Ae: 



4Gi:Ko{meg/r 



(55) 



Thus taking Lc ^ H and Cg being a fraction Cgo of the local 
scale height, we obtain Aest 4\n{l/{megoh))/{7vQ). On 
account of the logarithmic factor, the condition Aest > 1 
suggests that edge disturbances which lead to instabilities 
can be present for Q significantly larger than unity, as is 
confirmed numerically. However, a surface density threshold 
must be exceeded. Thus such modes associated with vorten- 
sity maxima do not occur in a non self- gravitating disc. Fi- 
nally we remark that although we applied a model assump- 
tion to obtain an integral equation with symmetric kernel 
in the above analysis, which enables the existence of solu- 
tions to be shown, the demonstration of a surface density 
threshold does not depend on this. The application of the 
Cauchy-Schwartz inequality to (|50)) may be carried out in a 
similar manner for the integral equation (|43|) leading to sim- 
ilar conclusions. The fundamental quantity is the vortensity 
profile, (|54)) indicates that if it is not sufficiently peaked, 
there can be no mode. Similarly, mode existence becomes 
less favoured for increasing softening and/or increasing m. 

At this point we remark that an analysis of the type 
discussed above does not work for vortensity minima. An 
equation similar to ()49|) could be derived but in this case 
the corresponding kernel IZ would be negative, implying in- 
consistent negative eigenvalues A. This situation results from 
the surface density perturbation leading to a gravitational 
potential perturbation with the wrong sign to satisfy the 
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condition (|4ip . Thus we do not expect the type of edge dis- 
turbance considered in this section to be associated with 
vortensity minima. The fact that edge modes associated with 
vortensity maxima require a threshold value of to be 
exceeded separates them from vortex forming modes which 
are modes associated with vortensity minima and require 
not to be above a threshold in order to be effective. For 
most disc models considered in this paper, with minimum 
Q ^ 2, vortex modes do not occur. Disc s with minimum 
Q > 2 are considered in lLin Papaloizoul (|201ll ). 

5.1.1 Implications for disc-planet systems 

For the locally isothermal disc models adopted in this paper, 
with a small aspect-ratio h = 0.05, edge modes could appear 
more easily as the disc is cold. This is because if the disc 
thickness sets the length scale of the edge, equations ()52|) . 
(|54p and ()55|) indicate that, for fixed Q, lowering sound- 
speed and hence the disc aspect ratio increases Aest- 

We can apply these equations to gaps opened by a Sat- 
urn mass planet, which is considered in linear and nonlinear 
numerical calculations later on. Without instabilities, these 
gaps deepen with time. There is also vortensity mixing in 
the co-orbital region as fluid elements pass through shocks 
and repeatedly execute horseshoe turns close to the planet. 
In a fixed orbit this reduces the gap surface density and the 
magnitude of the edge vortensity peaks. In this way the con- 
ditions for edge modes become less favourable with time. In 
this case edge modes are expected to develop early on dur- 
ing gap formation, if at all. However, this effect may be less 
pronounced if the orbit evolves because the planet migrates. 

For larger planetary masses such as Jupiter, a deeper 
gap will be opened but stronger shocks are also induced, 
which may lead to stronger vortensity peaks. Thus, in view 
of potentially competing effects, the conditions for edge 
modes as a function of planet migration and planetary mass 
must be investigated numerically and this will be undertaken 
in future studies. 



5.2 Launching of spiral density waves 

Although localised at the gap edge, the edge disturbance 
perturbs the bulk of the disc through its gravitational 
potential exciting density waves. This is expected to be 
through torques exerted at Lindblad resonances (Goldreich 
& Tremaine 1979). When the disc is exterior (interior) to 
the edge the wave excitation can be viewed as occurring at 
the the outer (inner) Lindblad resonance respectively. These 
resonances occur where a/m = —Q^K/m. Here the nega- 
tive (positive) sign alternative applies to the outer (inner) 
Lindblad resonance respectively. The perturbing potential is 
given by ()42p and from now on this potential is given the 

symbol ^edge- 

The total conserved angular momentum flux associated 
with the launched waves, when the y are assumed to prop- 
agate out of the system is given by iGoldreich Tremaind 
(|l979h as 

TT^mrH dredge 2mQa^edge ^ 



Ft 



dr 



+ 



(56) 



where (5 — 2n(^n' — mQ'). These waves carry positive an- 
gular momentum outwards when excited by an inner edge. 



or equivalently negative angular momentum inwards when 
excited by an outer edge. Thus they will destabilise nega- 
tive angular momentum disturbances at an inner disc edge 
or positive angular momentum disturbances at an outer disc 
edge that cause their emission. They will lead to an angu- 
lar momentum flux at the boundary where r — Vh given by 
(27r) [r{FGr + Fav)] |r, = Ft 



5.3 Spiral density waves and the growth of edge 
modes 

We now investigate the effects of wave losses at the non- 
sharp boundary as a perturbation. To do this we relax the 
assumption that the surface density tapers to zero at the 
boundary where r — Vh. Instead we adopt a small value 
there such that self-gravity is not important for the density 
waves. Thus we retain the domain V for evaluating the grav- 
itational potential to be (n,ro). Then, for real frequencies, 
all the terms in equation (|36|) apart from the term inversely 
proportional to a and the boundary term associated with 
the advective angular flux at r = are real. The imagi- 
nary part of the latter term was shown to be proportional 
to the wave angular momentum flux. We assume that the 
mode is close to marginal stability (subscript 'ms' below) 
with corotation at a vortensity maximum where a — Gr and 
write equation (|36)) in the form T{S^ a) = J^{S^ Gr + ^cr) = 0. 
Assuming ^cr is small, we then expand to first order in 5a ^ 
obtaining 

dT_ 
da 



5a + T{S, ar) = {n 

Differentiating (|36)) we find 

m\S\^ df] 



+ iD^)6a + J='{S,ar) = 0. (57) 



D 



and 



rfa'^ 



dr 



dr 



dS_ 

dr 



2mQa 



dr 



(58) 



(59) 



Setting 6a = 6ar — i'j, where 7 is the growth rate and taking 
the imaginary part of ([57)) . noting that the imaginary part 
of the first two terms on the right hand side of equation 
(|36)) contribute to the imaginary part of T giving this as 
Im(J^) = -=b(2/m) [r{FGr + Fav)] \r, = T(l/(m7r))FT with 
Ft given above, we obtain 

7 = -±-^+<5a.^. (60) 
Similarly the real part of (|57)) gives 

Drdar = -7 A - Re(J^), (61) 

where Re(J^) is the real part of T. We suppose that the 
waves emitted by the edge disturbance result in 7 7^ 0, but 
that the back reaction does not change the location of the co- 
rotation point from that of the marginally stable mode (as 
indicated by numerical results) . Thus co-rotation remains at 
the original vortensity maximum, implying that conditions 
adjust so that 5a r — 0. We then have: 

Ft 



mirDr 



(62) 
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where the upper (lower) sign apphes to an inner (outer) edge. 
In this case, instabihty occurs when wave emission causes 
negative (positive) angular momentum to be transferred to 
a negative (positive) angular momentum edge disturbance 
located at an inner (outer) sharp edge. 

According to equation (|62)) . to obtain instability (7 > 0) 
for an inner edge disturbance we need Dr < (since 
Ft > 0). For a disturbance concentrated at corotation near 
an inner sharp edge, with the disc lying beyond, the first 
term on the right hand side of (|58|) dominates. As corota- 
tion is at a vortensity maximum, we expect the contribution 
to the first integral of Dr from the region just beyond coro- 
tation, where dr/fdr < 0, to be negative and the contribution 
from the region just interior to corotation, where drj/dr > 0, 
to be positive. Because the region exterior to an inner edge 
has higher surface density than that interior to it, and the in- 
tegrand for the first term in Dr is proportional to S, we may 
expect the region exterior to co-rotation dominates the con- 
tribution to Dr, making it negative and therefore unstable. 
Indeed, we find Dr < for the numerical fiducial case pre- 
sented in section 16.11 In this case there is instability due to 
the reaction of the emitted outwardly propagating positive 
angular momentum wave on the negative angular momen- 
tum inner gap edge disturbance. A corresponding discussion 
applies to a disc with an outer sharp edge. 

We comment that the above considerations depend on 
the excitation of waves at a Lindblad resonance that were 
transported with a conserved action or angular momentum 
flux towards a boundary where they were lost. However, our 
linear calculations and simulations described below indicate 
lack of sensitivity to such a boundary condition. This can 
be understood if it is emphasised that the significant issue 
is that after emission, the wave action density should not 
return to the edge disturbance responsible for it. This is 
possible for example if the waves become trapped in a cavity 
in the wider disc in which case a radiative boundary would 
not be needed. 



6 LINEAR CALCULATIONS 

In this section, we present numerical solutions to the linear 
normal mode problem where the background self- gravitating 
disc contains a gap, presumed to have been opened by a 
planet. The basic state is obtained from hydrodynamic sim- 
ulations by azimuthally averaging the surface density and 
azimuthal velocity component fields to obtain one dimen- 
sional profiles. The basic state is thus axisymmetric. The 
radial velocity is assumed to be zero. 

We adopt a local isothermal equation of state P = 
C5(r)S with sound-speed Cs = hy^ GM^/r and h — 0.05 
for consistency with nonlinear simulations. The softening 
prescription used is that presented in 321 As for the ana- 
lytic discussion, the gravitational potential due to the planet 
and viscosity are neglected. The linearised equations follow 
from equations (|30|) - (|33|) together with equations and 
(|35p . However, rather than solve in terms of the quantity 
S' — E^Cs/SH-^^, which was convenient for analytic discus- 
sion, we find it more convenient here to work in terms of the 
relative surface density perturbation, W — for which 

a single governing equation may be written down explicitly. 
In terms of this, the velocity perturbations u'r^ can be 



written in the form 



1 

D 



^^=D 



2dW d^' 

^c^ I c,— \ — 

dr dr 



, dW d^' 



+ 



dr 



dr 



ma 
r 



(63) 
(64) 



The gravitational potential perturbation can be ex- 
pressed in terms of W through equations ()34p and ()35|) and 
we note that = r~^d{r^Q^) /dr. The derivatives d^' /dr 
and d^(^' /dr^ can be computed by replacing Km{r^r') with 
dKm/dr and Km/ dr^ ^ respectively. Inserting (|63)) and 
(|64p into the linearised continuity equation yields the gov- 
erning equation for W 



dr 



+ 



^2 dW ^ <W_ 
dr dr 



'r^ 
D 

2m d fi:Q 



+ 



2m 
a 



~D~ 



del 
dr 



■rT. 



W 



rD 



\W + ^') = C[W)^0. (65) 



Note that the the term in dc^/dr has been kept for con- 
sistency with simulations used to setup the basic state. We 
checked that this term has negligible effect on the results 
obtained below, by solving the linear problem without this 
term. This is because Cs varies on a global scale, whereas the 
edge disturbance is associated with strong local gradients. 

To solve equation (|65|) we discretised it on an equally 
spaced grid applied to the domain r = [1,10]. In prac- 
tice this employed Nr = 1025 grid points. Equation (|65|) 
together with the applied boundary conditions (see be- 
low) is then converted into a matrix equation of the form 
X^^-L £ij(cr)Wj = 0, where the matrix (Cij(a)) gives the 
discretised linear operator, which is a function of the fre- 
quency cr, and Wj is an approximation to W at the jih grid 
point. This problem cannot be solved for any value of a. 
This has to be consistent with the condition that the deter- 
minant of the system of linear equations be zero. Thus a is 
an eigenvalue although it is not an eigenvalue of (Cij) in the 
conventional sense. 

We proceed by taking (Cij) to be a function of a. For 
a specified value of a we solve the usual eigenvalue problem 
Cij{a)Wj = /j{a)Wi for an eigenvalue, /i, which may also 
be considered to be a function of a. The Newton-Raphson 
method is then used to solve the equation /i(cr) = 0. The val- 
ues of a so obtained are the required eigenvalues associated 
with the physical normal modes of the system. 

Because simulations suggest the important disturbances 
are associated with the outer gap edge, it is reasonable to 
search for values of a such that co-rotation lies near the outer 
gap edge. We checked that the reciprocal of the final matrix 
condition number is small (at the level of machine preci- 
sion) in order for results to be accepted. For simplicity, the 
boundary condition dW/dr = was applied at r = = 1 
and r = To = 10. As discussed in section 14. 6[ modes are 
found to be driven by the back reaction of emitted den- 
sity waves on a disturbance located at an outer gap edge, a 
process expected to be insensitive to boundary conditions. 
Indeed, we find that our results are insensitive as to whether 
the above or other boundary conditions are used (see section 
[631) . 



Edge modes 13 



2.5 



2.0 



1 .5 



1 .0 
0.5 
0.0 

4 5 6 7 

r 

Figure 4. Gap profile produced by a Saturn-mass planet in 
the Qo = 1-5 disc with viscosity v = 10~^. The surface density 
Ti and vortensity 77 are scaled by their values at r = 5.5. The 
relative deviation of the angular velocity from the Keplerian value 
is also plotted. Vertical lines indicate co-rotation radii rc (where 
Re (a) = 0). For the m = 2 mode Tc = 5.5, close to a vortensity 
maximum. When self- gravity is neglected Tc = 5.8, close to a 
vortensity minimum. 

6.1 A fiducial case 

In order to provide a fiducial case, we study the disc model 
with a gap, for which Qo = 1.5, that is adopted in JTl How- 
ever, the simulations described in JTlfrom which the model 
was extracted had u = 10~^, whereas our linear calculations 
described below are for inviscid discs. 

The basic state gap profile is illustrated in Fig. [H where 
the surface density S, the relative deviation of the angular 
velocity from the Keplerian value Q/Qk — 1 and the vorten- 
sity 7] are plotted. As expected from our analytic discussion 
of modes near to marginal stability, local extrema in 77 in the 
vicinity of gap edges are closely associated with instability. 
This is manifest through mode corotation points being very 
close to them. 

The magnitude of relative surface density perturba- 
tion, for self-gravitating (SG) and non- self- gravitating 
(NSG) responses are shown in Fig. \5\ For NSG cases we 
set = 0. The eigenfrequencies for SG (NSG) modes are 
given by -a = 0.1587 + i0.4515 x 10"^ {-a = 0.1458 + 
i0.2041 X 10~^). These correspond to co-rotation points at 
Tc = 5.4626, 5.7959 for the SG and NSG modes, respectively. 
The SG growth rate corresponds to ~ 3 times the local or- 
bital period. Thus although ^/(tr ^ 0.03, is relatively small, 
the instability grows on a dynamical timescale. 

The co-rotation points of SG and NSG modes are very 
close to local maximum and minimum of 77, respectively (Fig. 
ID). The SG mode grows twice as fast as the NSG mode, 
consistent with the observation made when comparing edge 
modes and vortex modes in 331 where the former became 
non-linear sooner. 

The SG and NSG eigenfunctions are similar around co- 
rotation (r G [5,6]), although the SG mode has a larger 
width and is shifted slightly to the left (Fig|5|). As expected 
from the discussion in section [5T2l the SG mode has signif- 
icant wave-like region interior to the inner Lindblad reso- 
nance (r ^ 3.4) and exterior to the outer Lindblad reso- 
nance at (r ^ 7.2). By contrast, the NSG mode has negligi- 
ble amplitude outside [5, 7] compared to that at co-rotation, 
whereas the SG amplitude for r G [8, 10] can be up to 56% 
of the peak amplitude near co-rotation. Increasing m in the 
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Figure 5. m = 2 unstable modes found from linear stability anal- 
ysis for the fiducial model, with self-gravity (solid) and without 
self-gravity (dotted). \W\ has been scaled by its value at corota- 
tion. 

NSG calculation increases the amplitude in the wave regions, 
but even for m = 6, we found the waves in [8, 10] for the NSG 
case have an amplitude of about 33% of that at co-rotation, 
a smaller value than that pertaining to the SG case with 
m = 2. 

The behaviour in the wider disc, away from co-rotation, 
shows that the for low m, the NSG mode is a localised vortex 
mode whereas the SG mode corresponds to an edge mode 
with global spirals. Noting that maxima in the vortensity 
and Toomre Q nearly coincide, it makes sense that the SG 
mode is global because away from co-rotation, the back- 
ground Toomre Q is decreasing, which makes it easier to 
excite density waves, because the evanescent zones between 
corotation and the Lindblad resonances that are expected 
from WKBJ theory, narrow accordingly. 

Comparing the SG and NSG modes shown in Fig. [5l 
we see that including self-gravity in the linear response en- 
ables additional waves in the disc at low m. Fig. [6] shows 
the gravitational potential perturbation for the SG mode. A 
comparison with in Fig.[5]shows that the surface density 
perturbation around Vc has an associated potential pertur- 
bation that varies on a more global scale. The peak in 
about Tc is confined to r G [4.6, 6.2] whereas that for |$^| is 
confined to r G [3.2,7.5], overlapping Lindblad resonances 
(see Fig. [2]) in the latter case. Thus is less localised around 
co-rotation than \W\. 

Rapid oscillations seen in \W\ for r > 7 are not ob- 
served for \^'\ in the same region. Furthermore, $^(r > 7) 
is at most :^ 20% of the co-rotation amplitude, which is a 
smaller ratio than that for \W\. This leads to the notion that 
the disturbance at the outer gap edge is driving the distur- 
bance in the outer disc as in our analytical discussion given 
above. In effect, the disturbance at co-rotation acts like an 
external perturber (e.g. a planet) to drive density waves in 
the outer disc, through its gravitational field. Clearly, this 
is only possible in a self- gravitating disc. 

6.2 Energy balance of edge and wavelike 
disturbances 

The analysis in JS] describes the edge mode as as being as- 
sociated with a coupling between disturbances associated 
with vortensity extrema near an edge and the smooth re- 
gions interior or exterior to the edge. We apply this idea 



y /y 














14 Lin and Papaloizou 




Figure 6. Gravitational potential perturbation (= in the 
plot), scaled by its value at r = 5.5 for the m = 2 mode with 
self- gravity. 



to the reference case (with self-gravity) by considering the 
region r > Vp. Specifically, we focus on r G [5, 10] because 
hydrodynamic simulations indicate the spiral arms are more 
prominent in the outer disc (33]and 3Zl)- 

Multiplying the governing equation ( I65p by S"" , inte- 
grating over [ri, r2] and then taking the real part, we obtain 



Re / gdr = Re / (^corot + g^ave)dr, 

J ri J r\ 

where 



^wave 



2m 



rE f2dW ^ d^ 
D \ ^ dr dr 



dr 

2m /EQX dd 



D J dr 



m E |cv/|2 
rD 



(66) 

(67) 
(68) 

(69) 



g has dimensions of energy per unit length and, for a normal 
mode, its real part is four times the energy per unit length 
associated with pressure perturbations (the term oc ) and 
gravitational potential perturbations (the term oc $^*). As 
the factor of four is immaterial to the discussion we simply 
call the integral of this quantity over the region concerned 
the thermal-gravitational energy (TGE). We remark that 
when self- gravity dominates, the TGE is negative and when 
pressure dominates, it is positive. 

When the integral is performed, one sees that the TGE 
is balanced by various terms on the RHS of equation ()66|) . 
For simplicity we split the terms on the RHS into just two 
parts that are integrals of ^corot and ^wave over the region 
of interest. This is of course not a unique procedure. The 
vortensity term ^corot has been isolated because it contains 
the potential co-rotation singularity which can be amplified 
by large vortensity gradients at the gap edge. The rest of the 
RHS is collected into ^wave- Note that the term in dc^/dr 
is included m ^wave 5 

despite being proportional to 1/cr. This 
is motivated by trying to keep ^corot as close to the vorten- 
sity term identified in the analytical formalism as possible. 
However, we have considered attributing the dc^/dr term to 
^corot or neglecting it altogether. In both cases it made neg- 
ligible difference compared to the splitting adopted above. 



Again, this is because varies on a much larger scale than 
vortensity gradients. 

It is important to note that while the real part of the 
combination ^corot + ^wave gives rise to the TGE, we can not 
interpret Re(^corot) as an energy density of the co-rotation 
region. It contains a term contributing to the TGE that 
is proportional to the vortensity gradient (see below) and 
potentially associated with a corotation singularity. Similar 
arguments apply to ^wave which, in the strictly isothermal 
case, can be seen to be associated with density waves (see 
section above) . We use this splitting to show that for the 
modes of interest, the vortensity term contributes most to 
the TGE. 

Numerically however, quantities defined above are in- 
convenient because of the vanishing of D for neutral modes 
at Lindblad resonances. To circumvent this we work with 
D^g, D^gcoTotj and D^g^a^e- We call these modified energy 
densities. This change does not disrupt our purpose because 
the most important balance turns out to be between the 
terms involving Re{D^ g) and Re(D^pcorot), both focused 
near co-rotation, where D ^ . In the region near Lindblad 
resonances and beyond, Re(^) and Re(^corot) are small. Thus 
the incorporation of the factor does not influence conclu- 
sions about the TGE balance. The modified energy densities 
are plotted in Fig. [71 for the m = 2 mode in the fiducial case. 
The curves share essential features with those obtained us- 
ing the original definitions without the additional factor of 

(these are presented and discussed in Appendix lA)) . 

Fig. [71 shows that Re{D^ g) is negative around co- 
rotation which is located near the outer gap edge, r G [5, 6]). 
Accordingly Re(^) < and must be dominated by the grav- 
itational energy contribution since the pressure contribution 
is positive definite. This supports the interpretation of the 
disturbance as an edge mode where self-gravity is impor- 
tant. If it were the vortex mode then Iie{D^ g) > near 
corotation, because in that case, self- gravity is unimportant 
compared to pressure. 

For r ^ 8.4, Re(D^p) becomes positive due to pressure 
effects and oscillates towards the outer boundary as the per- 
turbation becomes wave-like. Similarly, the pressure pertur- 
bation dominates towards the inner boundary (not shown). 
This signifies that self-gravity becomes unimportant rela- 
tive to pressure within these regions. This is consistent with 
the behaviour of the gravitational potential perturbations, 
which are largest near co-rotation. We checked explicitly 
that the gravitational energy contribution to the TGE is 
largely focused around the gap edge (r G [5,6]), even more 
so than the gravitational potential perturbation. 

Fig. [3 shows that Re(D^pcorot) and Re(D^pwave) have 
their largest amplitudes for r G [5,6]. Integrating over r G 
[5, 10], we find U < 0, and using a normalisation such that 



i 



U = Re D^gdr= -\U\, 



we find 



= Re 



r 
i: 



D'gcorotdr -0.822\U\ 



D g^ 



edr : 



This implies the TGE 
gravitationally-dominated 
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disturbance. 



and hence a 
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Figure 7. One-dimensional modified energy densities computed 
from the eigenfunctions W, for the fiducial case Qo = 1.5. 



conclude that if co-rotation occurs at a maximum value of 77 
then the vortensity term contributes negatively to the TGE. 
As stated above, when self- gravity dominates pressure, the 
TGE is negative, so it can then be mostly accounted for by 
the vortensity term. This is precisely the type of balance 
assumed for the neutral edge mode modelled in section [5Tl 
However, a consequence of the above discussion is that if 
co-rotation occurs at a minimum value of 77, then the vorten- 
sity term contributes positively to the TGE. This can only 
give the dominant contribution when the TGE is positive, 
which can only be the case when the pressure contribution 
dominates over that from self-gravity. In the limit of weak 
self-gravity, the TGE will be > and it can be balanced 
by a localised disturbance with co-rotation at a vortensity 
minimum as happens for the vortex modes. 



\Ucorot/U\ 0.8, suggesting the TGE is predominantly 
balanced by the vortensity term which from Fig. [71 is 
localised in the gap edge region, balances the gravitational 
energy of the mode. Because vortensity gradients are largest 
near the gap edge, gravitational energy is most negative 
here overtaking pressure, resulting in a negative TGE. 

6.2.1 Further analysis of the vortensity term 

The vortensity term pcorot that we defined above differs from 
that naturally identified in the alternative form of the gov- 
erning equation in 311 However, this difference is not signif- 
icant. The vortensity term can be further decomposed as 



^corot — ^corot,l ~h ^corot,2 

d (\ 



^corot,l — 



m 



^corot,2 



cr(l — zy^) dr \r\ 
2m dv 



\S' 



\S' 



where v — a I k. Let us now temporarily consider ^corot,i as 
the new vortensity term, so that ^corot ^corot,i and at- 
tribute ^corot,2 to the wave term so that ^wave ^wave + 

^corot,2. The new vortensity term is proportional to the 
vortensity gradient explicitly. With the new definitions we 
find 

|t/corot/t/| - 0.690. 

Thus the new co-rotation term alone accounts for ^ 70% 
of the (modified) TGE and giving almost the same result 
as the previous one. In addition as most of the contribution 
comes from around co-rotation, where \v\ <^ 1. Replacing 



the factor (1 



by unity has a negligible effect on the 



Pcorot,! 



results. However, making this replacement gives 

rn\S^d_ rV 
a dr \ri ^ 

which corresponds to the vortensity term used in the analy- 
sis given in 311 Whether this term contributes positively or 
negatively to the TGE depends on the sign of ^corot,l dr. 
We expect most of the contribution to this integral comes 
from co-rotation where a 0. Then sgn (^f^^ gcorot,idr^ = 

sgn( ^corot,i 1^ ). Evaluating this for a neutral mode with co- 
rotation at a vortensity extremum, we find ^corot,i|^^ = 
-r]~'^\S'\'^(d'^r]/dr'^)/Q'\ . Since quite generally, < 0, we 



6.3 Dependence on m and disc mass 

We have solved the linear eigenmode problem for discs with 
1.2 ^ Qo ^ 1.6. The basic state is set up from hydrodynamic 
simulations as described in 331 with v = 10~^. The local 
max((5) near the outer gap edge ranges between Q = 3.3 and 
Q = 4.4. Self-gravity is included in the response / 0) . 

FigOa) compares eigenfunctions , \ W\, for different m 
for the Qo = 1-2 disc. Increasing m increases the ampli- 
tude in the wave region relative to that around co-rotation 
(rc ^ 5.5), resulting in m = 6 being more global than m = 3. 
High m modes do not fit our description of edge modes in 
the region of interest (r > 5), because the disturbance at co- 
rotation becomes comparable to or even smaller than that 
in the wave region beyond the outer Lindblad resonance. It 
is then questionable to interpret the waves as a secondary 
phenomenon induced by the edge disturbance, even though 
the physics of the modes, as being entities driven by an un- 
stable interaction between edge and wave disturbances, may 
be more or less the same. Hence, we typically find what we 
describe as edge-dominated modes with m < 3 in simu- 
lations where the surface density perturbation is maximal 
near the gap edge. For fixed m = 2, hence in that regime, 
Fig. [8jb) shows that lowering Qo makes the co-rotation am- 
plitude larger relative to that in the wave region. This is 
expected because increasing the level of self-gravity means 
the necessary condition for edge disturbance is more easily 
satisfied (section ISTT]) . 

In Fig. [9lwe plot the locations of the corotation points 
and the growth rates for the unstable modes we found as a 
function of the azimuthal mode number, m and a range of 
disc masses (parametrised by Qo)- Note that changing Qo 
affects both the background state and the linear response, 
while m acts as a parameter of the linear response only. The 
plots in Fig. [9ja) show that the co-rotation radius tends to 
move outwards with increasing m and/or Qo. However, coro- 
tation is always located in the edge region where the vorten- 
sity is either decreasing or maximal. The largest growth rates 
are found for low m (< 3) edge-dominated modes. 

Referring to the basic state (Fig. (H), we see that in- 
creasing m and/or Qo, which weakens self-gravity, shifts co- 
rotation towards max(?7~^) (the vortensity minimum), which 
disfavours edge dominated modes. For sufficiently low mass 
or non- self- gravitating discs, we expect co-rotation associ- 
ated with vortensity minima and therefore vortex modes 
to dominate (e.g. Qo = 4 in 33]). Thus we only expect 
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Figure 8. The effect of azimuthal wave- number m (a) and the 
effect of disc mass, parametrised by Qo (b) on hnear edge modes. 



Figure 9. Co-rotation radii (a) and growth rates (b) as a function 
of azimuthal wave- number m, for a range of Qo values. 



edge dominated modes for low m and sufficiently large disc 
mass. An exception to the general trend above is that the 
Qo = 1.6, m = 6 mode, has a smaller co-rotation radius 
than the corresponding mode with Qo < 1.6. This may be 
a boundary effect because for this value of m and Qo the 
modes are distributed through the outer disc and require 
the implementation of accurate radiative boundary condi- 
tions, rather than the simplified boundary conditions actu- 
ally applied. However, such high m modes typically have 
growth rates smaller than those for lower m, accordingly we 
do not expect them to dominate, nor are they observed in 
the nonlinear simulations discussed later. 

Figinib) shows that growth rates increase as Qo is low- 
ered (increasing surface density scale). We found that de- 
creasing Qo results in deeper gaps, because disc material 
trapped in the vicinity of the planet adds to the planet po- 
tential, so that the effective planet mass is larger for smaller 
Qo. Steeper gaps are expected to be more unstable, there- 
fore there is a contribution to the increased growth rate from 
changes to the background profile as the disc mass is in- 
creased. 

We also expect the effect of lowering Qo, through the 
linear response, to de-stabilise edge dominated modes as 
they rely on self- gravity. However, the effect of self- gravity 
through the response is weaker for larger m values because 
the size of the Poisson kernel decreases. Hence, the increase 
of the growth rates with disc mass for m = 4 — 6 is likely 
more attributable to the effect of self-gravity on the ba- 
sic state. Growth rates for the most unstable mode ap- 
proximately doubles as the disc mass is increased by 25% 
(Q^ = 1.5 ^ Qo = 1.2). Notice for fixed Qo ^ 1.5, the 
m = 3 — 5 growth rates are similar but m > 3 are not edge 
modes. 

For Qo = 1.2, 1.3 the most unstable mode has m = 3, 
whereas for Qo = 1.5, the m = 2 mode has the largest 



growth rate. As discussed in ^ edge dominated modes re- 
quire sufficient disc mass and/or low m. For fixed m, they 
are stabilised with increasing Qo. Hence, for edge modes to 
exist with decreasing disc mass, they must shift to smaller 
m. The higher m modes are less affected by self- gravity. 
This may explain the double peak in growth rate plotted as 
a function of m for Qo = 1.6 as we move from Qo = 1.5 
(we remark that for Q^ = 1.5 the m = 3 growth rate is also 
slightly smaller than m — 2, 4). 

The integrals of the modified energy densities for each 
azimuthal wave-number are given for the fiducial case with 
Qo = 1.5 in Table [U Only the modes with m ^ 2 are edge 
dominated modes because for these sum of the integrals con- 
tributing to the (modified) TGE is negative, corresponding 
to a gravitationally dominated disturbance and over 50% 
of the energy is accounted by the vortensity edge term. The 
mode with m — 2 had the largest growth rate and corotation 
radius at r = 5.46. This is fully consistent with the nonlin- 
ear simulation of the fiducial case performed in section [3Jl 
There the dominant mode in the outer disc had m — 2 and 
corotation radius at r ^ 5.5. 

For m > 2, the total (modified) energy becomes pos- 
itive, which must be due to the pressure term {cl\i:'\/T) 
becoming dominant. The vortensity term alone cannot bal- 
ance this because it contributes negatively. For m ^ 4, the 
vortensity term contribution is small consistent with high m 
solutions being predominantly wave- like (Fig. [8]). However, 
note that inspite of this the vortensity term may still play 
some role in driving the modes through corotation torques. 
These trends have also been found for other models so they 
appear to be general. However, high m global modes are 
unimportant for the problems we consider, as they are not 
seen to develop in nonlinear simulations. 
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Table 1. Eigenfrequencies expressed in units of flk(r = 1), coro- 
tation radii and modified energy integrals for different azimuthal 
mode numbers, m, for the fiducial disc with Qo = 1-5. The mod- 
ified energy integrals are taken over the outer disc r £ [5, 10]. 



m -CTR 7 X 10^ rc -Ucorot/\tj\ -f/wave/|f/| 



1 


0.079 


0.270 


5.454 


0.65 


0.35 


2 


0.159 


0.452 


5.463 


0.82 


0.18 


3 


0.237 


0.404 


5.498 


0.34 


-1.24 


4 


0.313 


0.412 


5.533 


0.026 


-1.03 


5 


0.386 


0.366 


5.585 


0.0048 


-1.004 


6 


0.462 


0.242 


5.603 


0.0015 


-0.998 


7 


0.524 


0.223 


5.699 


0.0010 


-0.997 


8 


0.599 


0.189 


5.699 


0.00099 


-0.995 



6.4 Softening length 

Gravitational softening has to be used in a two-dimensional 
calculation. It prevents a singularity at r = in the Poisson 
kernel and approximately accounts for the disc's vertical di- 
mension but is associated with some uncertainty. Apart from 
the linear response, softening also has an effect through the 
gap profile set up by our nonlinear simulations. A fully self- 
consistent treatment requires a new simulation with each 
new softening considered. For reasons of numerical tract abil- 
ity though, we only performed experiments using two values 
of softening parameter e^o,e in simulations to setup the gap 
profile and range of values e^o,i used in the linear response. 

Note that in nonlinear simulations a single disc poten- 
tial softening parameter e^o is used. Our results are sum- 
marised in Table [2] The integrated total modified energy 
values indicate we have found dominated edge modes since 
the major contribution is due to the co-rotation/vortensity 
term. 

Comparing growth rates for fully self- consistent cases 
e^o,e = e^o,i = 0.3, 0.6 shows that increasing the softening 
length stabilises edge modes, which is expected since they 
are driven by self-gravity. In addition all growth rates for 
epO,e = 0.6 are smaller than those for e^ce = 0.3 indepen- 
dently of e^o,i- This indicates stabilisation of edge modes via 
the basic state when softening is increased. For fixed e^o,e, 
7 decreases as e^o,i increases. We could not find edge modes 
for e^o,i ^ 0.6 with e^o,e = 0.3, nor for e^o,! ^ 0.83 with 
epO,e = 0.6. An upper limited is expected from analytical 
considerations (see ^5.1|) because if self-gravity in the lin- 
ear response is too weak, the vortensity /co-rotation distur- 
bance cannot be maintained so edge dominated modes are 
suppressed. This is intuitive because the edge mode requires 
gravitational interaction between gap edge and the smooth 
disc. 



6.4-1 Convergence issues 

For fixed e^o,e = 0.3, the energy ratios in Table [2] indicate 
convergence as e^o,i 0. Perhaps surprisingly, \Ucovot/U\ 
decreases with softening. We found this is because ^wave, 
as defined by equation (|69)) , includes a term proportional to 
dP^' /dr^ which is the dominant contribution to |?7wave|, and 
its contribution mostly comes from the co-rotation region (as 
the potential perturbation is largest there) . As the strength 
of self-gravity is increased by decreasing softening, the edge 



Table 2. As for tabled but for linear calculations with different 
softening parameters, for the Qo = 1.5 fiducial case and azimuthal 
mode number m = 2. Subscripts '1' denote the softening param- 
eter used in the linear response, and 'e' denotes the softening 
parameter used in setting up the basic state. Note that although 
growth rates vary by about a factor of two, the location of the 
corotation radius does not vary significantly. 





£(70,6 


-CTR 


7 X 10^ 


-f/corot/|C/| 


-t/wave/|C/| 


0.03 


0.3 


0.159 


0.560 


0.66 


0.33 


0.05 


0.3 


0.159 


0.544 


0.67 


0.32 


0.1 


0.3 


0.159 


0.508 


0.69 


0.30 


0.2 


0.3 


0.159 


0.461 


0.74 


0.23 


0.3 


0.3 


0.159 


0.452 


0.82 


0.18 


0.5 


0.3 


0.158 


0.435 


0.87 


0.14 


0.3 


0.6 


0.159 


0.420 


0.91 


0.09 


0.5 


0.6 


0.158 


0.410 


0.95 


0.064 


0.6 


0.6 


0.158 


0.375 


0.99 


0.0017 


0.8 


0.6 


0.158 


0.276 


1.15 


-0.14 



disturbance is modified, but the subsequent effect of the 
d^(^' /dr^ term on the energies cannot be anticipated a priori. 

For 6^0,1 = 0.03, the vortensity term does not domi- 
nate the modified total energy as much. However, this case 
is in fact not self-consistent because the softening used to 
setup the gap profile is an order of magnitude larger than 
that used in linear calculations. The limit of zero softening 
is also physically irrelevant because the disc has finite thick- 
ness. For the self- consistent cases with reasonable softening 
values (0.3, 0.6), the vortensity term dominates the energy 
balance, so the analytic description developed above works 
reasonably well. 



6.5 Edge mode boundary issues 

The boundary condition W' = was applied for sim- 
plicity. Vortex modes in low mass or non-self-gravitating 
discs are localised and in sensitive to boundary conditions 
(|de Val-Borro et ahllioOTl ). The edge dominated modes rely 
on self-gravity and are therefore intrinsically global, bound- 
ary effects cannot be assumed unimportant a priori. 

We performed additional experiments with varying 
boundary conditions for the fiducial case [Qo — 1-5). 
These include re-locating the inner boundary to r = 
1.1, 2.5 or the outer boundary to r = 9.8 and approximat- 
ing/ext rapolating bound ary derivatives using interior grid 
points (|Adams et al.|[l989l In the last case, the linear ODE 
is applied at endpoints of the grid and therefore no bound- 
ary conditions imposed. For self- gravitating disc calculations 
giving edge dominated modes, these various conditions gave 
co-rotation radii varying by about 0.2% and growth rates 
varying by at most 10%. Overall the eigenfunctions W are 
similar, showing the essential features of the edge mode, in- 
cluding relative amplitude of outer disc wave region and the 
co-rotation region. We conclude that the existence of edge 
dominated modes is not too sensitive to boundary effects. 
Physically this is because edge modes are mainly driven by 
the local vortensity maximum or edge, which is an internal 
feature away from boundaries. 
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7 NONLINEAR HYDRODYNAMIC 
SIMULATIONS 

We present nonlinear hydro dynamic simulations of disc- 
planet models with 1.2 ^ Qo ^ 2, corresponding to disc 
masses 0.079M* ^ ^ 0.047M*. Most simulations use a 
viscosity of v = 10~^ or equivalent ly a = O(10~^), which 
has been typically adopted for protoplanetary discs. This 
value has also been foun d to suppress vortex instabilities 
(|de Val-Borro et al.1 lioOTl ). We adopt gravitational soften- 
ing parameter Cgo = 0.3. The planet is set on a fixed cir- 
cular orbit at Tp = 5 in order to focus on gap stability. We 
briefly explore the effects of varying viscosity and softening 
lengths later in this section and planetary migration in the 
next section. 



7.1 Numerical method 

The hydro dynamic equations are evolved with the FARGO 
code (jMassetllioOQl ). FA RGO is an explicit finite- difference 
code similar to ZEUS (|Stone k NormanI Il992l ) but cus- 
tomised for disc-planet interactions. It circumvents the time 
step limit imposed by the rotational velocity at the inner 
boundary by splitting the azimuthal velocity into mean and 
perturbed parts, azimuthal transport being performed on 
each separate ly. Self-gravity for FARGO was implemented 
and tested bv lBaruteau Massetl (|2008l V Two-dimensional 
self-gravity can be calcula ted using Fast Fourier Transform 
(|Binnev &: TremaindfToS?! ). This requires the radial domain 
to be logarithmically spaced and doubled in extent. The 
planet potential is introduced at 25Po and its gravitational 
potential ramped up over lOPo, where Pq is the Keplerian 
period at r = 5. 

The disc is divided into NrXN^ = 768 x 2304 grid points 
in radius and azimuth giving a resolution of Ar/H = 16.7. 
The grid cells are nearly square, with Ar/rA(f) =1.1. We im- 
pos e open bound aries at ri and non-reflecting boundaries at 
To (|Godonlll996l ). The latter has also been used in the self- 



gravitating disc-planet simulations of I Zhan g et al.' {"2008') . 
As argued above, because edge modes are physically driven 
by an internal edge, we do not expect boundary conditions 
to significantly affect whether or not they exist. Indeed, sim- 
ulations with open o uter boundaries, or damping boundaries 
(|de Val-Borrdl2006l ) all show development of edge modes. 



7.2 Overview 

We first present an overview of the effect of edge dominated 
modes on disc-planet systems. Fig. [10] shows the relative 
surface density perturbation (AS/E) for Qo = 1.2, 1.5 and 
2.0. The profile formed in these cases has local max((5) = 
3.3, 4.2, 5.3 near the outer gap edge and the average Q for 
r G [6, To] is 1.6, 2.0, 2.6. Although the edge mode is asso- 
ciated with the local vortensity maximum (located close to 
max(Q) for gap profiles), it requires coupling to the external 
smooth disc via self- gravity. Hence, edge modes only develop 
if Q in these smooth regions is sufficiently small, otherwise 
global disturbances cannot be constructed. 

The case Qo = 2.0 gives a standard result for gap- 
opening planets, typically requiring planetary masses com- 
parable or above that of Saturn. For that mass a partial 
gap of depth ^ 50% is formed, a steady state attained and 



remains stable to the end of the simulation. This serves as 
a stable case for comparison with more massive discs. This 
simulation was repeated with u = 10~^, in which case we 
identified both vortex and edge modes. The standard vis- 
cosity u = 10~^ thus suppress both types of instabilities for 
Qo = 2. Although we adopted v = 10~^ to avoid complica- 
tions from vortex modes, we note that they are stabilised for 
sufficiently massive discs anyway (e.g. Qo = 1.5, = 10~^ 
do not show vortices, see section [3|). 

Gap edges become unstable for Qm ^ 1.5 when u — 
10~^. As implied by linear analysis, the m — 2 edge mode 
dominates when Qo = 1.5 and it saturates by t = IOOPq. It is 
more prominent in the outer disc, and over-densities deform 
the gap edge into an eccentric ring. The under-density in 
the gap also becomes more non-axisymmetric, being deeper 
where the gap is wider. The inner disc remains fairly circu- 
lar, though non-axisymmetric disturbance (m = 3) can be 
identified. Note the global nature of the edge modes in the 
outer disc: spiral arms can be traced back to disturbances at 
the gap edge rather than the planet. Consider the spiral dis- 
turbance at the outer gap edge just upstream of the planet 
where AS/S is locally maximum. Moving along this spiral 
outwards, AS/S reaches a local minimum at (x, y) — (6, —4) 
before increasing again. This is suggestive of the outer disc 
spirals being induced by the edge disturbance. 

When the effect of disc self-gravity is increased to the 
Qo — 1.2 case, there is significant disruption to the outer 
gap edge, it is no longer clearly identifiable. The m = 1 edge 
mode becomes dominant in the outer disc, although overall 
it still appears to carry an m = 2 disturbance because of 
the perturbation due to the planet. The inner disc becomes 
visibly eccentric with an m = 2 edge mode disturbance. The 
shocks due to edge modes can comparable strength to those 
induced by the planet. 



7.3 Development of edge modes for Qo = 1.5 

We examine the fiducial case with Qo = 1.5. Fig. [TTI shows 
the development of the edge instability. The final dominance 
of m = 2 is consistent with linear calculations. The planet 
opens a well-defined gap by t 40Po, or ^ 5Po after the 
planet potential has been fully ramped up. At t 

46Po, two over-dense blobs, associated with the edge mode 
can be identified at the outer gap edge. At this time the 
gap has AS/S ~ —0.3, implying edge modes can develop 
during gap formation, since the steady state stable gap in 
Qo = 2 reaches AS/S ~ —0.6. Edge modes are global and 
are associated with long trailing spiral waves with density 
perturbation amplitudes that are not small compared to that 
at the gap edge. This is clearly seen at t 50Po when spiral 
shocks have already developed. We estimate that the edge 
modes have a growth rate consistent with predictions from 
linear theory and become non-linear within the time frame 
46Po < t < 50Po. 

The edge perturbations penetrate the outer gap edge 
and trail an angle similar to the planetary wake. Notice the 
m = 3 disturbance near r = 4 in the snapshots taken at 
t = 50Po and t = 56Po. Two of the three over-densities, 
seen as a function of azimuth, correlate to the m — 2 mode 
in the outer disc, while the third over-density adjacent to 
the planet correlates with the outer planetary wake. These 



Edge modes 19 



99.95 orbits 



99.95 orbits 



99.95 orbits 



Figure 10. Surface density perturbations (relative to t = 0) at t ~ lOOPo as a function the minimum value of Toomre Q in the disc, 
Qo. From left to right: Qo = 2.0, 1.5, 1.2. 



features result from self-gravity and are unrelated to vortex 
modes. 

The presence of both planetary wakes (with pattern 
speed ^^fc(5)) and edge modes (with pattern speed :^ 
nfc(5.5)) of comparable amplitude enable direct interaction 
between them. Passage of edge mode spirals through the 
planetary wake may disrupt the former, explaining some of 
the apparently split-spirals in the plots. At t ^ 64Po Fig. [TT] 
indicates that a spiral density wave feature coincides with 
the (outer) planetary wake. The enhanced outer planetary 
wake implies an increased negative torque exerted on the 
planet at this time. This effect occurs during the overlap of 
positive density perturbations. Assume that at a fixed ra- 
dius the edge mode spiral has azimuthal thickness nh where 
n is a dimensionless number. The time taken for this spiral 
to cross the planetary wake is AT = nh/\Qk{^-^) — ^k{^)\ — 
O.OGnPo. This is ^ Pq for n = 0(1). We expect such cross- 
ing spiral waves to induce associated oscillations in the disc- 
planet torque. 

The non-linear evolution of edge modes in disc-planet 
systems can give complicated surface density fields. The gap 
can become highly deformed. Fig. [TT] shows that large voids 
develop to compensate for the over-densities in spiral arms, 
leading to azimuthal gaps (t = 64Po). At t = 72Po, the 
gap width ahead of the planet is narrowed by the spiral 
disturbance at the outer edge trying to connect to that at 
the inner edge. However, this gap-closing effect is opposed 
by the planetary torques which tend to open it. 

A quasi-steady state is reached by t ^ 76Po showing an 
eccentric gap. For t = 99.5Po additional contour lines are 
plotted to indicate two local surface density maxima along 
the edge mode spiral adjacent to the planet. This spiral is 
disjointed around r ^ 7.5 where on traversing an arm, there 
is a minimum in AS/E. Taking co-rotation of the spiral at 
Tc = 5.5 gives the outer Lindblad resonance of an m = 2 
disturbance at r = 7.2, assuming a Keplerian disc. This is 
within the disjointed region of the spiral arm, which sup- 
ports our interpretation that the disturbance at the edge is 
driving activity in the outer disc by perturbing it gravita- 
tionally and launching waves at outer Lindblad resonances. 



7.3.1 Eccentric gaps 

Several snapshots of the Qo = 1.5 run taken during the sec- 
ond half of the simulation are shown in Fig. 1121 which display 
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Figure 12. Precession of an eccentric gap in the presence of an 
m = 2 edge mode for a disc with Qo = 1.5. 



the precession of an eccentric outer gap edge. Deformation 
of a circular gap into an eccentric one requires an m = 1 
disturbance. By inspection of the form of the relative sur- 
face density perturbation for Qo = 1.5 (see Fig. [TO]) . we see 
that edge modes are associated with an eccentric outer gap 
edge. However, the inner gap edge remains fairly circular. 

In their non-self-gravitating disc-planet simulations, 
iKlev DirksenI (|2006l ) found eccentric discs can result for 
planetary masses Mp > 0.003M* fixed on circular orbit. Our 
simulations show that edge modes in a self- gravitating disc 
can provide an additional perturbation to deform the gap 
into an eccentric shape for lower mass planets. 



7.4 Evolution of the gap structure for Qo = 1.2 

With increasing disc mass, evolution of the gap profile away 
from its original form takes place. Fig. [13] illustrates the 
emergence of the m = 3 edge mode, expected from linear 
calculations, and subsequent evolution of the gap profile. 
By t ^ 42Po the m = 3 edge mode has already become 
non-linear, whereas for Qo = 1.5 the m — 2 mode only 
begins to emerge at t ^ 44Po, being consistent with linear 
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Figure 11. Evolution of the surface density perturbation (relative to t = ) for the Qo = 1.5 case from t ~ 40Po to t ~ IOOPq- 



calculations that the former has almost twice the growth 
rate of the latter. In Fig. [13] there is a m = 3 spiral at 
r < Tp in phase with the mode for r > Vp. The gap narrows 
where they almost touch. These are part of a single mode. 
By t = 44Po, shocks have formed and extend continuously 
across the gap. The self- gravity of the disc has overcome the 
planet's gravity which is responsible for gap-opening. Edge 
mode spirals can be more prominent than planetary wakes. 

The evolution of the form of the gap is shown in the 
one-dimensional averaged profiles plotted in Fig. [13] The 
quasi-steady profile before the onset of the edge mode in- 
stability is manifest at t ~ 40Po. By t — 44Po, the original 
bump at r = 6 has diminished. This bump originates from 
the planet expelling material as it opens a gap but is subse- 
quently 'undone' by the edge mode as it grows and attempts 
to connect across the gap. This gap filling effect is possible 
as non axisymmetric perturbation of the inner disc may be 
induced via the outer disc self-gravity. Notice also in Fig. 
[T3] at t = 44Po the increased surface density at r ^ 7.5 
implying radial re-distribution of material. By t — 120Po, 
the inner bump (r = 4) is similar to that produced by a 
standard gap-opening planet. By contrast, the outer bump 
cannot be maintained. The gap widens and there is an over- 
all increase in surface density for r > 6.5. Note that there 
is no additional diminishing of the outer edge bump from 
t 44Po 120Po. 




Figure 13. The disc model with Qo = 1-2: emergence of the 
m = 3 edge mode (top) and evolution of the gap profile (bottom). 
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Figure 14. The disc model with Qo = 1.2: evolution of the 
azimuthally averaged vortensity 77 profile in the region of the gap. 



Fig. [Ml illustrates the evolution of the vortensity profile 
within the gap. The development of edge modes temporarily 
reduces the amplitude of the vortensity maxima set up by 
the perturbation of the planet. This is particularly notice- 
able when the profiles at t = 40Po and t = 44Po are com- 
pared. However, vortensity is generated through material 
passing through shocks induced by the gravitational pertur- 
bation of the planet ( Lin & Papaloizou 2010). This provides 
a vortensity source that enables the amplitudes of vortensity 
maxima to increase again, as can be seen from the profile 
at t = 6OP0. The maxima remain roughly symmetric about 
the planet's location at rp ± 2rh. Note a general increase 
in the co-orbital vortensity level caused by fluid elements 
repeatedly passing through shocks. 

Although the m = 3 mode emerges first and is predicted 
to be most unstable from linear analysis, it does not persist. 
Linear theory can only predict the initially favoured mode. 
In the nonlinear regime, Fourier analysis shows that m — 2 
becomes dominant for 6OP0 ^ ^ ^ 145Po . We suspect this be 
due to finite viscosity acting over this time-scale. Larger m 
means shorter radial wavelengths, so over given a time-scale, 
diffusion and hence stabilisation by viscosity is more effec- 
tive than smaller m. The evolution is also complicated by 
the planetary wake which may gravitationally interact with 
edge modes. For Qo = 1.2 we found that m = 1 becomes 
dominant at t > 145Po when there appears to be a coupling 
between the outer planetary wake and the edge mode. 

7.5 Mass transport 

The similarity between the effects induced by an edge mode 
and an external perturber (planet) is further illustrated by 
the induced mass transport. Just as when a planet opens a 
gap, the development of an edge disturbance results in a ra- 
dial re-distribution of mass. We express the angular momen- 
tum transport in terms of an a viscosity parameter, defined 
through a — AurAu^/cl, where A here denotes deviation 
from the azimuthal average. 

Angular momentum is also transmitted through gravi- 
tational torques. However, in practice we found this to be a 
small effect compared to transport due to Reynolds stresses. 
The discs here are not gravito-turbulent. The azimuthally 
averaged a parameters (a) are shown in Fig. [151 The sta- 
ble disc with Qo = 2 has non-axisymmetric features entirely 
induced by the planetary perturbation and (a) is localised 
about the planet's orbital radius giving a two- straw feature 



about rp. This is typical of disc-planet interactions, and it 
provides a useful case for comparison with more massive 
discs to see the effect of edge modes and their spatial de- 
pendence. 

Increasing self- gravity by adopting Qo — 1.5, we see 
that angular momentum transport is enhanced around rp 
and for r < 4, signifying the global nature of edge modes. 
Towards the inner boundary, a ^ 10~^ being twice as large 
as for the case with Qo = 2. The Qo = 1.5 case also has 
enhanced wave-like behaviour for r > 7 as compared to Qo = 
2.0. 

The Qo = 1.2 case is more dramatic, with (a) reaching a 
factor of 2-3 times larger than the imposed physical viscosity 
around the planet's location (z///i^ = 4 x 10~^). This case 
displays a three- straw feature, with an additional trough at 
about r = 5.5, i.e. close to edge mode co-rotation, as if 
another planet were placed there. Although the Qo = 1.5 
case also develops edge modes, there is no such additional 
trough at co-rotation. We suspect this may because the Qo = 
1.2 case is more unstable, developing an m = 3 mode ( which 
produces 3 additional effective coorbital planets placed at 
the gap edge), whereas the Qo = 1.5 case develops a less 
prominent m = 2 mode. However, compared to Qo — 1.5, 
Qo = 1.2 has no enhanced transport away from the gap 
region. 

Fig. [151 also shows the evolution of the disc mass as a 
function of time, Mdit). Viscous evolution leads to mass loss 
(Qo = 2.0) and there is clearly enhanced mass loss if edge 
modes develop (Qo = 1.2, 1.5). For t ^ 50Po the curves are 
identical, this interval includes gap formation and the emer- 
gence of edge modes. When edge modes become non-linear, 
mass loss is enhanced (t :^ 50Po), but there is no signifi- 
cant difference between Qo = 1.2, 1.5, which is consistent 
with the previous {a) plots near the inner boundary. The 
Qo = 1.2 curve is non-monotonic, though still decreasing 
overall, possibly because of boundary effects, and mass loss 
is enhanced once more around t — 125Po. As Qo is lowered, 
mass loss becomes less smooth, showing the dynamical na- 
ture of edge modes. 

7.6 Additional effects 

We briefly discuss effects of softening and viscosity on edge 
modes. To analyse mode amplitudes, the surface density 
field is Fourier transformed in azimuth giving am{r) as the 
(complex) amplitude of the m^^ mode. We then integrate 
over the outer disc r > 5 to get Cm — f (imdr and consider 
A — Cm/ Co. 

7. 6. 1 Softening 

In linear theory we found that edge modes cannot be con- 
structed for gravitational softening parameters that are too 
large. We have simulated the Qo = 1.5 disc with e^o = 
0.6,0.8, 1.0.0 

Fig. [161 shows the evolution of the running-time aver- 
aged m — 2 amplitudes. The evolution for all cases is very 

^ It should be remarked that increasing e^o makes the Poisson 
kernel in simulations increasingly non-symmetric, so the earlier 
analytical discussion is less applicable. 
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Figure 15. Top: azimuthally averaged Reynolds stresses asso- 
ciated with edge modes when they become non-linear (Qo = 
1.2, 1.5), the Qo = 2 case is shown as a stable case with no edge 
modes for comparison. Bottom: evolution of the total disc mass 
for the three cases, scaled by initial mass. 

similar for t ^ 40Po- Unstable modes develop thereafter, 
with increasing growth rates and saturation levels as Cgo is 
lowered. Cgo = 1 displays either no growth or a much re- 
duced growth rate. The gap for Qo = 1.5, Cgo = 1.0 reaches 
a steady state similar to the stable case with Qo = 2.0, Cgo = 
0.3. This is consistent with the fact that as softening weak- 
ens, edge disturbances can no longer perturb the remainder 
of the disc via self-gravity, the residual non axisymmetric 
structure being due to the planetary perturbation. 

Interestingly, the decrease in saturation level in going 
from Cgo = 0.3 0.6 is smaller than going from Cgo = 0.6 ^ 
0.8, which is in turn smaller than Cgo = 0.8 1.0, despite a 
larger relative increase in softening in one pair than the next. 
This is suggestive of convergence for the nonlinear saturation 
amplitude as ^0. This can be expected from convergence 
in linear theory f ^6.4.ip . However, at Cgo = 0.3 there are 
only 5 cells per softening length. To probe even smaller Cgo 
requires prohibitively high resolution simulations. Further- 
more, since Cg approximately accounts for the disc's non-zero 
vertical extent, very small softening lengths are not physi- 
cally relevant to explore. 

7.6.2 Viscosity 

An important difference between edge dominated modes and 
the already well-studied vortex modes in planetary gaps is 
the effect of viscosity on them. The standard viscosity ly = 
10~^ prevents vortex mode development, but we have seen 
in 33] that such viscosity values still allow the m — 2 edge 
modes to develop. Vortices are localised disturbances and 



-1 .5 










-2.0 










-2,5 




/ / 
/ •" ' --'^ 


6^0=0-3 

£go=0-6 


















-3-0 










25 


50 75 

i/orblts 


1 00 


i; 



Figure 16. Running-time average of the surface density m = 2 
Fourier amplitude, integrated over r G [5, 10] and scaled by the 
axisymmetric component, for the Qo = 1.5 disc, as a function of 
the gravitational softening length. 
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Figure 17. Non-axisymmetric disc modes for the Qo = 1.5 disc 
as a function of applied viscosity v. The Fourier amplitude has 
been integrated over r G [5, 10] and scaled by the axisymmetric 
amplitude and its running-time average plotted. Thick lines indi- 
cate the m = 2 mode and is plotted for all u. Thin lines indicate 
the m = 3 mode and is plotted for u = 10~^ — 10~^ only. 



more easily smeared out than global spirals in a given time 
interval. Hence, vortex growth is inhibited more easily by 
viscosity than spirals. 

We repeated the Qo = 1.5 runs with a range of viscosi- 
ties. Results are shown in Fig.[T71 Generally, lowering viscos- 
ity increases the amplitudes of the non-axisymmetric modes. 
For 1/ ^ 10~^, the m = 3 mode emerges first (note that this 
is not in conflict with our linear calculations which used 
u — 10~^ in its basic state), rather than m — 2. This is be- 
cause lowering viscosity allows a sharper vortensity peak to 
develop in the background model, and thus has the same ef- 
fect as increasing the disc mass. The latter can enable higher 
m edge modes. 

Unlike vortex modes, even with twice the standard vis- 
cosity (v — 2 X 10~^) the edge mode develops and grows. It 
is only is suppressed when ^ 5 x 10~^. However, this is be- 
cause no vortensity maxima could be setup at the gap edges 
due to vortensity diffusion in the co-orbital region making 
an almost uniform vortensity distribution in the gap. Since 
the necessary condition for edge modes is not achieved, no 
linear modes exist. This differs from the nature of the sup- 
pression of vortex modes, because in that case vortensity 
extrema that are stable can still be setup. 
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8 APPLICATIONS TO DISC-PLANET 
SYSTEMS 

In order to focus on the issue of the stabihty of the dip/gap 
in the surface density profile induced by a giant planet, 
the planet was held on a fixed orbit. However, torques ex- 
erted by the non axisymmetric disc on the planet will in- 
duce migration that could occur on a short time scale (eg. 
IPeplihski et aLlfioOSl V Such torques will be significantly af- 
fected by the presence of edge modes. Accordingly, we now 
discuss the torques exerted by a disc, in which edge modes 
are excited, on the planet. 

8.1 Disc-planet torques 

The presence of large-scale edge mode spirals of comparable 
amplitude as the planetary wake will significantly modify the 
torque on the planet, which in a stable disc is the origin of 
disc-planet torques. Here, we measure the disc-planet torque 
but still keep the planet on fixed orbit. 

Fig. [TST a) shows the evolution of the torque per unit 
length for the Qo = 1.5 fiducial case. The torque profile 
at t = 40Po, before the edge modes have developed signifi- 
cantly, shows the outer torque is larger in magnitude than 
inner torque, implying inward migration. This is a typical 
result for disc-planet interactions and serves as a reference. 
At t = 60Po and t = lOOPo, edge modes mostly modify the 
torque contributions exterior to the planet, though some dis- 
turbances are seen in the interior disc (r — < — 5r^, rn 
being the Hill radius). The original outer torque at +r^, is 
reduced in magnitude as material re-distributed to concen- 
trate around the co-rotation radius Tc of the edge mode, (rc 
is ^ 2rh away from the planet). 

Edge modes can both enhance or reduce disc-planet 
torques associated with planetary wakes. Consider r — G 
[3, 5]rh in Fig. llSf a). Comparing the situation at t = 60Po 
to that at t = 40Po , the torque contribution from this region 
is seen to be more negative. This is because an edge mode 
spiral overlaps the planetary wake and therefore contributes 
an additional negative torque on the planet. However, at 
t = lOOPo, an edge mode spiral is just upstream of the 
planet, exerting a positive torque on the planet, hence the 
torque contribution from r — Vp G [3, 5]r/i becomes positive. 

Differential rotation between edge modes and the planet 
produces oscillatory torques, shown in Fig. [TST b). Unlike 
typical disc-planet interactions, which produce inward mi- 
gration, the total instantaneous torques in the presence of 
edge modes can be of either sign. For a disturbance with 
m fold symmetry the time interval between encounters with 
the planet is AT = (27r/m)/^Q, where 50. is the difference 
in angular velocity of the planet and the disturbance pat- 
tern. Approximating 50 — \Ok{rp) — Qfc(rc)| with Vc — 5.46 
obtained from linear theory, we obtain AT = 4.04Po. In- 
deed, the oscillation period in Fig. [18] is :^ 4Po. Both inner 
and outer torques oscillate, but since the edge mode is more 
prominent in the outer disc, oscillations in the outer torque 
are larger, particularly after mode saturation (t > 70Po). 

We compare time-averaged torques in Fig. [TSl^ c). In all 
three cases, on average a negative torque acts on the planet 
and its magnitude largest at t = 40Po. Up to this point, the 
torque becomes more negative as the disc mass increases, as 
expected. However, development of edge modes makes the 



averaged torques more positive, and eventually reverses the 
trend with Qo. While the disc with Qo — 1.5 also attains 
steady torque values as in Qo = 2, the disc with Qo = 1.2 
has significant oscillations even after time- averaging and is 
remains non-steady at the end of the run. 

Before the instability sets in, torques arise from wakes 
induced by the planet, and the outer (negative) torque 
is dominant. Edge modes concentrate material into spiral 
arms, leaving voids in between. Therefore, except when spi- 
ral arms cross the planet's azimuth, the surface density in 
the planetary wake is reduced because it resides in the void 
in between edge mode spirals. Hence the torque magnitude 
is reduced. If the reduction in surface density due to edge 
mode voids is greater than the increase in surface density 
scale produced by decreasing Qo, the presence of edge modes 
will, on average make the torque acting on the planet more 
positive. 

8.2 Outward scattering by spiral arms 

If edge modes develop, planetary migration may be affected 
by them. Their association with gap edges inevitably affects 
co-orbital flows. While a detailed numerical study of migra- 
tion is deferred to future work, we highlight an interesting 
effect found when edge modes are present. This is outward 
migration induced through scattering by spiral arms. 

We restarted some of the simulations described above 
at t = 50Po allowing the planet to move in response to the 
gravitational forces due to the disc. The equation of motion 
of the planet is integrated using a fifth order Runge-Kutta 
integratojfl. Fig. [19] shows the orbital radius of the planet 
in discs with edge modes present. No obvious trend with 
Qo is shown. This is because of oscillatory torques due to 
edge modes and the partial gap associated with a Saturn 
mass planet, which is asso ciated with type HI migration 
(|Masset PaDaloizoull2003h . The initial direction of migra- 
tion can be inwards or outwards depending on the relative 
positioning of the edge mode spirals with respect to the 
planet at the time the planet was first allowed to move. 

Outward migration is seen for Q^ = 1.3 and Qo = 1.5. 
For Qo = 1.5, the planet migrates by Ar = 2, or 8.6 times 
its initial Hill radius, within only 4Po. This is essentially the 
result of a scattering event. Repeating this run with tapering 
applied to contributions to disc torques from within O.Gr/^ 
of the planet, we found outward scattering still occurs, but 
limited to Ar = 1 and the planet remains at :^ 6 until 
t = 128Po. Hence, while the subsequent inward migration 
seen for Qo = 1.5 may be associated with conditions close 
to the planet, the initial outward scattering is due to an 
exterior edge mode spiral. 

For Qo = 1.3 the planet is scattered to Vp :^ 8. It is 
interesting to note that the subsequent inward migration 
for Qo = 1.5 and Qo = 1.3 stalls at r ^ 6, i.e. the original 
outer gap edge. The planet remains there for sufficient time 
for both gap and edge mode development and for Qo = 

1.3 a second episode of outward scattering occurs (it also 
occurs for Qo = 1.5 to a small extent around t — 70Po). 

3 iLin Papaloizoul (|2Q10h found the same integrator to be ade- 
quate for studying migration induced through scattering by vor- 
tices 
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Figure 18. (a) Disc torque per unit radius acting on the planet. The length scale = (Mp/S^^^rp is the Hill radius, (b) Time 
dependent evolution of the disc planet torques, the contribution from the inner disc (dotted curve), the outer disc (dashed curve) and 
the total torque (solid curve), (c) The running time average of the total torque acting on the planet as a function of Qo- 
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Figure 19. Orbital radius evolution in discs with edge mode 
disturbances. The instantaneous orbital radius rp is shown as a 
function of time. 



Interaction with edge modes can affect the disc well beyond 
the original co-orbital region of the planet. In the case of 
outward migration, it may promote gravitational activity in 
the outer disc. However, boundary effects may be important 
after significant outward scattering. 

The spiral arm-planet interaction for Qo — 1.5 above 
is detailed in Fig. 1201 A spiral arm approaches the planet 
from the upstream direction (t = 50.9Po) and exerts positive 
torque on the planet, increasing Vp. The gap is asymmetric in 
the azimuthal direction with the surface density upstream of 
the planet being larger than that downstream of the planet 
(t = 51.5Po). As the spiral arm passes through the plan- 
etary wake (t — 51.5Po — 52.1Po) the gap surface density 
just ahead of the planet builds up as an edge mode is setup 
across the gap. Notice the fluid blob at r = 5, ^p — Lpp — O.Stt. 
This signifies material executing horseshoe turns ahead of 
the planet. At the later time t — 53.1Po the planet has 
exited the original gap, in which the average surface den- 
sity is now higher than before the scattering. Material that 
was originally outside the gap loses angular momentum and 
moves into the original gap. This material is not necessarily 
that composing the spiral. 

We remark that migration thr ough scattering by spir al 
arms differs from vortex scattering (|Lin PapaloizoulbOlO ). 
Vortices form about vortensity minima, which lie further 
from the planet than vortensity maxima. We can assume 
vortensity maxima are stable in the case of vortex formation. 
This means the ring of vortensity maxima must be disrupted 
in order for vortices to flow across the planet's orbital radius 
for direct interaction. Edge modes themselves correspond to 
a disruption of vortensity maxima, hence direct interaction 
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Figure 20. Interaction between an edge mode spiral and the 
planet (green dot), causing the latter to be scattered outwards, 
log E is shown. 



is less hindered than in the vortex case. Furthermore, a vor- 
tex is a material volume of fluid. Vortex-planet scattering 
results in an orbital radius change of both objects. A spiral 
pattern is generally not a material volume of fluid. It can be 
seen in Fig. 1201 that the local max(S) in the spiral remains 
at approximately the same radius before and after the in- 
teraction. The spiral as a whole does not move inwards. In 
this case, the spiral first increases the planet's orbital radius, 
thereby encouraging it to interact with the outer gap edge 
and scatter that material inwards. 



9 SUMMARY AND DISCUSSION 

We have studied instabilities associated with a surface den- 
sity gap opened by a giant planet embedded in a massive 
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protoplanetary disc. Vortex producing instabilities, associ- 
ated with local vortensity minima, have previously been 
found to occur in weakly or non self- gravitating discs. How- 
ever, edge modes that are associated with local vorten- 
sity maxima and require sufficiently strong self-gravity were 
found in the more massive discs that we have focused on in 
this paper. These have very different properties in that they 
are global rather than local and are associated with large 
scale spiral arms as well as being less affected by viscous 
diffusion. 

For our disc-planet models with fixed planet mass Mp = 
3 X 10~^M* and fixed disc aspect-ratio h = 0.05 we found 
edge modes to develop for gap profiles with an average 
Toomre Q < 2 exterior to the planet's orbital radius. In 
the unperturbed smooth disc, this corresponds to Q < 2.6 
at the planet's location and Q < 1.5 at the outer boundary, 
or equivalently a disc mass Md > 0.06M* . For fixed M^, non- 
linear simulations show edge modes develop readily for uni- 
form kinematic viscosity z/<2xl0~^ and self- gravitational 
softening Cg < 0.8H. 

A theoretical description of edge modes was developed. 
The edge mode is interpreted as a disturbance associated 
with an interior vortensity maximum that requires self- 
gravity to be sustained, and which further perturbs the re- 
mainder of the disc through its self-gravity causing the ex- 
citation of density waves. Linear calculations confirm this 
picture and also show the expected strengthening of the in- 
stability when self-gravity is increased via the disc mass and 
the expected weakening of the instability through increasing 
gravitational softening. 

Hydrodynamic simulations of disc-planet interac- 
tions were performed . We confirm earlier suggestions by 
iMeschiari &; LaughlinI (|2008l ) , who found gravitational insta- 
bilities associated with their prescribed gap profiles without 
a planet. Indeed, we found that edge modes can develop 
for gaps self-consistently opened by introducing a planet 
into the disc. However, our models showed their develop- 
ment during gap formation of a Saturn-mass planet, when 
the gap only consists of a 20 — 30% deficit in surface den- 
sity relative to the unperturbed disc. This is much less 
shallow than the Jovian-ma ss planetary gaps considered in 
IMeschiari LaughlinI (|2008l l , which are typically associated 
with a 90% deficit. We found edge modes to exist for disc 
models extending to a distance twice as large as the plan- 
ets orbit with masse s Md 0.06M^, again this is less mas- 
sive than required in IMeschiari LaughlinI (|2008i ). but more 
massive than those typically used in modelling protoplane- 
tary discs. 

We remark that in IMeschiari LaughlinI ' s disc model, 
the local Q maximum occurs at the gap centre, presumably 
where the planet would lie. However, gap opening by the 
planet may disrupt the disturbance associated with this ex- 
tremum, which is required to induce perturbations in the 
smooth disc. It is then unclear if the edge mode could be 
setup. Furthermore, since corotation lies at their gap cen- 
tre, there would be no relative motion between spirals and 
the planet of the kind seen in our simulations. The case of 
Jovian-mass planets will be explored in a future work. We 
note in that case, significant di sruption of the disc, lead- 
ing to fragmentati on, may occur (|Armitage &: Hansenll 19991 : 
iLufkin et al1l2004l l. 

Our simulations show edge modes with m = 2 and 



m — 3. They have surface density maxima localised near 
the outer gap edge, but the disturbance they produce ex- 
tends throughout the entire disc, consistent with analytical 
expectation and linear calculations. One important differ- 
ence between edge modes and the vortex forming modes in 
non self- gravitating discs is that the latter require low viscos- 
ity. The typical viscosity values adopted for protoplanetary 
discs, ly — 10~^, will suppress vortex formation but not edge 
modes associated with self-gravity. 

We have also considered the effect of edge modes on 
planetary migration. They produce large oscillations in the 
disc torques acting on the planet, which can be positive or 
negative. The presence of edge modes reverses the trend of 
the time averaged disc-on-planet torque as a function of disc 
mass: the torque being more positive as disc mass increases. 
Direct interaction between the planet and a spiral arm asso- 
ciated with an edge mode is possible. These should be more 
prominent in the disc section with the lower Q. In practise 
this corresponds to r > r^, so the planet tends to interact 
with the outer spirals, which results in a scattering of the 
planet outwards. 

9.1 Outstanding issues 

This work has been motivated by the wish to obtain a better 
understanding of disc-planet interactions in massive discs. 
In particular a proper discussion of phenomena such as type 
HI migration requires incorporation of self-gravity which we 
have undertaken here. We have studied model discs which 
are unstable through the development of edge modes when a 
dip/gap is produced by an embedded giant planet but which 
would be stable without the planet. 

Our first calculations show that the disc torques acting 
on the planet in the presence of edge modes change sign and 
are highly time-dependent. Thus migration is unpredictable. 
A statistical approach may be ultimately required to assess 
the likelihood of reduced inward migration in practice. Fur- 
thermore, development of edge modes during gap formation 
indicates their importance for planets before they reach a 
Jovian mass. This leads to the possibility of type HI migra- 
tion, which is self-sustained and can be in either direction. 
Issues such as what is the appropriate softening length to 
use in a two dimensional simulation and the treatment of 
material inside the Hill radius, will be important for more 
thorough understanding of planetary migration in the pres- 
ence of large-scale spiral arms in massive discs. Resolving 
these requires improved numerical modelling which will be 
presented in a future work. 
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APPENDIX A: ENERGY DENSITIES IN 
LINEAR THEORY 

In ^we defined and discussed the thermal- gravitational en- 
ergy (TGE) and the different contributions to it. There, a 
factor was applied to the TGE per unit length and to 
each contributing term in order to overcome numerical dif- 
ficulties associated with Lindblad resonances. Por complete- 
ness, we provide here a discussion of the behaviour of the 
TGE and the various contributions to it considered without 
the additional factor . 

Pig. lAll shows very smilier features to the correspond- 
ing curves discussed in ^6.21 Re(^) is negative around co- 
rotation, again signifying a self-gravity driven edge mode. 
Beyond r :^ 6.4, Re(p) becomes positive due to pressure. 
The most extreme peak at r ^ 5.6 coincides with the 
background vortensity edge (Pig.|4]). The negative contribu- 
tion from the co-rotation term is larger in magnitude than 
that from the positive wave term, resulting in Re(p) < 0. 
The co-rotation radius rc < 5.6 making the shifted fre- 
quency ^^(5.6) = m[Q{5.6) — n(rc)] < 0. Also at r = 5.6, 
d(ri~^)/dr > 0, resulting in Re(pcorot) < at this point. 
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Figure Al. The thermal and gravitational energy (TGE) per 
unit length (solid) computed from eigensolutions for the fiducial 
case Qo = 1-5 from linear theory. The contributions from the 
vortensity term (dotted line) and other terms (dashed line), de- 
fined by Eq. [67l are also shown. The relatively small bumps near 
r = 7 are numerical. However as they are the same magnitude 
for both contributions, incorporating them does not affect con- 
clusions concerning the overall energy balance in the system. 

Re(pwave) and Re(pcorot) have relatively small spurious 
bumps around r = 7.2 that are associated with the outer 
Lindblad resonance and are numerical. The eigenfunctions 
W, are well-defined without singularities there, but eval- 
uation of pcorot and pwave luvolvcs divlslou by D, which 
can amplify numerical errors at Lindblad resonances where 
D ^ 0. Integrating Re(p) over [5, 10], we find [/ < 0, the 
TGE is negative, which means gravitational energy domi- 
nates over the pressure contribution for r ^ 5. Integrating 
the contributions to the TGE separately, we find 

t/corot = Re / pcorotdr -0.94|[/|, 

PwaveC^r ^ -0.0771 [/|. 

The integration range includes the OLR and thus the spu- 
rious bumps in pwave and pcorot- However, we still find that 
U approximately equals t/corot + tAvave- The correct energy 
balance is still maintained despite being subject to possible 
numerical error due to the diverging factor 1/D. The above 
means that, aside from the spurious bumps, the energy den- 
sity values for r G [5, 10] may still be used to interpret energy 
balance. We find \Ucorot/U\ ^ 0.9, so as before, the TGE is 
predominantly accounted for by the vortensity term. 



